ioam / topographica

A general-purpose neural simulator focusing on topographic maps.
topographica.org
BSD 3-Clause "New" or "Revised" License
53 stars 32 forks source link

How to change the geometry of v1 area ? #656

Closed windform2015 closed 8 years ago

windform2015 commented 8 years ago

Hi, I would like to ask, how to change the geometric shape of v1 area ? In GCAL model, the default area of v1 is a square. Now I am kind of exploring what is the dynamics of v1 activity like if the v1 has another type of geometry, e.g. a rectangle instead of a square. I suppose it requires some changes in the parameters in ''param.Number()" and I am in the process of looking for that. Could you give a small hint about this ? Thank you.

jbednar commented 8 years ago

There's a mask parameter for the V1 sheet that can be set to whatever shape you like.

jbednar commented 8 years ago

Oh, but if you just want a rectangle, you don't have to use a mask; just double the width of each sheet (Retina, LGN, and V1). A mask will be wasteful of memory in that case. The easiest way to do that is to edit examples/gcal_oo_or.ty to change the nominal_bounds_template for each of those sheets, replacing radius=X with points=((-2X,-X),(2X,X)).

You'll also need to change anything that's calculated from the area, which is then no longer valid (since you're manually doubling the area). That means the total_num_inputs will need to be doubled (since there's twice the retinal area to cover), and the lbound and ubound on the x parameter of the inputs (in inputs=...) will need to be doubled, so that the input patterns are chosen from the full area of the retina.

So it's easy enough to do, but you have to check that you do everything correctly. You could of course automate this by adding an area_aspect_ratio parameter and then calculate everything from that, but doing it manually as above is probably less work.

windform2015 commented 8 years ago

Thanks for so detailed explanation !

windform2015 commented 7 years ago

Hi, up to now I tested by doubling both areas of each sheet and the connection fields between them, but it still doesn't work. Then I realized probably the density of each sheet should be also changed, so I set : for example,

topo.sim['Retina']=sheet.GeneratorSheet(xdensity=48,ydensity=24, input_generator=retina_inputs, period=1.0, phase=0.05, nominalbounds=sheet.BoundingBox(points=((-2(p.area/2.0+0.25+0.375+0.5),-(p.area/2.0+0.25+0.375+0.5)),(2_(p.area/2.0+0.25+0.375+0.5),p.area/2.0+0.25+0.375+0.5))))

there are errors shown.

I am not sure the direction is correct. Could you help have a look ?

jbednar commented 7 years ago

If you want to change the area, you shouldn't also change the connection field areas and the density; neither of those have anything to do with the area. See my 2004 Neuroinformatics paper for the detailed reasoning, but if you scale everything, including the connection fields, you'll just end up right back where you started, having simply increased the density (neurons simulated per unit area) in one direction rather than the area (extent of the cortical surface that you are modelling).

So, do not change the CF extents, and do not change the densities, if you want to change the area. Just look in gcal_oo_or.ty for anything that is calculated from p.area -- these are the only things you should be changing, because these are the ones that area affected when the area changes. Personally, I'd leave p.area set to its current value (1.0 by default), but then just change how it's used as appropriate to reflect that you will now have a rectangular sheet instead of a square one. E.g. it's used in the calculation of total_num_inputs-- if you're doubling the width of the sheet, multiply this equation by 2.0 to make it correct now. The next thing is inputs -- you'll want to change the area that the inputs are selected from, making the width twice as long as before. Then comes the retina's size -- you'll double the width, and leave the other dimension the same. Same for the rest of the sheets -- each of them uses p.area, and you'll need to change each one to reflect the fact that you now want a rectangle.

To actually do this, you have three options:

  1. You can go through and replace "p.area" with "2.0*p.area" everywhere in gcal_oo_or.ty when it's being used for a width, while not changing the specifications that have to do with the length.
  2. Or if you wanted to do it really cleanly, you'd make a new parameter sheet_aspect_ratio and replace "p.area" with "p.aspect_ratio*p.area" wherever you want a width, and use just "p.area" as before whenever you want a length.
  3. Or you could replace "p.area" with two parameters "p.length_scale" and "p.width_scale", and replace every "p.area" in the file with one of those two parameters, depending on whether it's being used for a width or a length.

So far, we've only needed one parameter "p.area", since we use square sheets, but if you want non-square sheets you effectively need two parameters to control sheet sizing, one for width and one for length. Hope that makes sense!

And I hope it works, which it should. If it doesn't, post your diffs (the file before and after your changes), along with a detailed error message.

windform2015 commented 7 years ago

Hi , thanks for your detailed answer ! I tested that as your instruction and I chose size(p.area)=1.5 and cortex_density=49, for the matrix of orientation preference I got, it is 73×73 instead of 146 ×73. Is this because I did something wrong or the default data area for shown is chopped to size(p.area)^2 ? I checked python files, e.g. topo.base.simulation, topo.base.sheetviews, topo.base.sheet and so on, but haven't figured this out. Sorry to disturb you again !

jbednar commented 7 years ago

You don't need to change cortex_density, but 49 should be fine for it.

I don't know what you mean by size(p.area); what is "size"? And can you clarify which of the three options above you are using (1, 2, or 3)? What you are describing doesn't sound like any of them. There is no cropping or "default data area" that I think would be relevant here.

windform2015 commented 7 years ago

By size(p.area), I mean area=param.Number(default=1.5,bounds=(0,None), inclusive_bounds=(False,True)). Based on the the first option you suggested, I changed it as: e.g. In the original script, topo.sim['Retina']=sheet.GeneratorSheet(nominal_density=p.retina_density, input_generator=retina_inputs, period=1.0, phase=0.05, nominal_bounds=sheet.BoundingBox(radius=p.area/2.0+0.25+0.375+0.5)) By changing nominal_bounds, it becomes topo.sim['Retina']=sheet.GeneratorSheet(nominal_density=p.retina_density, input_generator=retina_inputs, period=1.0, phase=0.05, nominal_bounds=sheet.BoundingBox(points=((-(2.0*p.area/2.0+0.25+0.375+0.5),-(p.area/2.0+0.25+0.375+0.5)),((2.0*p.area/2.0+0.25+0.375+0.5),p.area/2.0+0.25+0.375+0.5))))

In summary, I used 2.0*p.area to replace "p.area" everywhere in gcal_oo_or.ty, and meanwhile, changed nominal_bounds=sheet.BoundingBox(radius=X) to nominal_bounds=sheet.BoundingBox(points=((-2X, -X),(2X,X))).

Then from the simulation result, for the matrix of orientation preference I got, it is 73×73 instead of 146 ×73. However, the area of v1 has become a rectangle area with the width which is a double of length. I am not sure, the data I got is right ? Please have a look.

Also I have another puzzle, for example, for the Retinal sheet, could the nominal_bounds be changed to nominal_bounds=sheet.BoundingBox(points=((-2.0*(p.area/2.0+0.25+0.375+0.5),-(p.area/2.0+0.25+0.375+0.5)),(2.0*(p.area/2.0+0.25+0.375+0.5),p.area/2.0+0.25+0.375+0.5))) ?

Thanks for your attention !

jbednar commented 7 years ago

I'm not sure precisely what you did from that description, but it sounds correct, except that you can't replace p.area everywhere, only where it is being used to scale the width. Changing from radius to points is correct, and seems like it's written properly. But beyond that I'm not sure what's happening, so please post your actual .ty file so that I can take a look.

windform2015 commented 7 years ago

Okay, sorry for not showing you very clearly.

windform2015 commented 7 years ago

try: import external except: pass

import sys, topo

if not (sys.version_info[0] == 2 and sys.version_info[1] == 7): print "Warning: Topographica requires Python 2.7, and will not currently work with any other version." sys.exit()

from topo.misc.commandline import process_argv

import argparse parser = argparse.ArgumentParser() parser.add_argument('-inputseed')

parser.add_argument('-cortex_density')

parser.add_argument('-aff_learning')

parser.add_argument('-retina_density')

parser.add_argument('-lgn_density')

args = parser.parse_args()

print parser.parse_args()

from math import pi import os

import numpy import param import imagen from imagen import deprecated from topo import learningfn,numbergen,transferfn,pattern,projection,responsefn,sheet

import topo.learningfn.optimized import topo.learningfn.projfn import topo.transferfn.optimized import topo.pattern.random import topo.pattern.image import topo.responsefn.optimized import topo.sheet.lissom import topo.sheet.optimized

import topo.transferfn.misc from topo.base.arrayutil import DivideWithConstant from topo.transferfn.misc import HomeostaticResponse from topo.misc.commandline import global_params as p from holoviews.core import SheetCoordinateSystem

p.add(

input_seed = param.Integer(default=int(args.inputseed),doc="""
     The random seed to set the input patterns"""),

dataset=param.ObjectSelector(default='Gaussian',objects=
    ['Gaussian','Nature', 'VGR'],doc="""
    Set of input patterns to use::
      :'Gaussian': Two-dimensional Gaussians
      :'Nature':   Shouval's 1999 monochrome 256x256 images.
      :'VGR':       Simulated vertical goggle rearing
                   (anisotropically blurred Shouval)"""),

num_inputs=param.Integer(default=4,bounds=(1,None),doc="""
    How many input patterns to present per unit area at each
    iteration, when using discrete patterns (e.g. Gaussians)."""),

aff_strength=param.Number(default=1.5,bounds=(0.0,None),doc="""
    Overall strength of the afferent projection to V1."""),

exc_strength=param.Number(default=1.7,bounds=(0.0,None),doc="""
    Overall strength of the lateral excitatory projection to V1."""),

inh_strength=param.Number(default=1.4,bounds=(0.0,None),doc="""
    Overall strength of the lateral inhibitory projection to V1."""),

aff_lr=param.Number(default=0.1,bounds=(0.0,None),doc="""
    Learning rate for the afferent projection to V1."""),

exc_lr=param.Number(default=0.0,bounds=(0.0,None),doc="""
    Learning rate for the lateral excitatory projection to V1."""),

inh_lr=param.Number(default=0.3,bounds=(0.0,None),doc="""
    Learning rate for the lateral inhibitory projection to V1."""),

area=param.Number(default=1.5,bounds=(0,None),
    inclusive_bounds=(False,True),doc="""
    Linear size of cortical area to simulate.
    2.0 gives a 2.0x2.0 Sheet area in V1."""),

retina_density=param.Number(default=24,bounds=(0,None),
    inclusive_bounds=(False,True),doc="""
    The nominal_density to use for the retina."""),

lgn_density=param.Number(default=24,bounds=(0,None),
    inclusive_bounds=(False,True),doc="""
    The nominal_density to use for the LGN."""),

cortex_density=param.Number(default=49,bounds=(0,None),
    inclusive_bounds=(False,True),doc="""
    The nominal_density to use for V1."""),

retinal_waves=param.Integer(default=0,bounds=(0,None),doc="""
    How many retinal wave (noisy disk) presentations before the
    chosen dataset is displayed. Zero for no noisy disk
    presentations otherwise 6000 has been typically used.."""),

percent_noise = param.Number(default=50, doc="""
   The percentage of the noise strength to the disk strength for noisy disks"""),

contrast=param.Number(default=100, bounds=(0,100),
    inclusive_bounds=(True,True),doc="""
    Brightness of the input patterns as a contrast (percent)."""),

gain_control = param.Boolean(default=True, doc="""
    Whether or not gain-control (divisive lateral inhibitory) is to be
    applied in the LGN"""),

homeostasis = param.Boolean(default=True, doc="""
    Whether or not the homeostatic adaption should be applied in V1"""),

t_init = param.Number(default=0.20, doc="""
    The initial V1 threshold value. This value is static in the L and GCL models
    and adaptive in the AL and GCAL models.""")

)

retinalbounds = sheet.BoundingBox(points=((-2(p.area/2.0+0.25+0.375+0.5),-(p.area/2.0+0.25+0.375+0.5)),(2_(p.area/2.0+0.25+0.375+0.5),(p.area/2.0+0.25+0.375+0.5))))

contrast_scale = p.contrast / 100.0 if p.dataset=="Gaussian": total_num_inputs=int(p.num_inputs_p.area__2) inputs=[pattern.Gaussian(x=numbergen.UniformRandom(lbound=-(2.0_p.area/2.0+0.25), ubound= (2.0*p.area/2.0+0.25),seed=p.input_seed+12+i), y=numbergen.UniformRandom(lbound=-(p.area/2.0+0.25), ubound= (p.area/2.0+0.25),seed=p.input_seed+35+i), orientation=numbergen.UniformRandom(lbound=-pi,ubound=pi,seed=p.input_seed+21+i), size=0.088388, aspect_ratio=4.66667,scale=contrast_scale) for i in xrange(total_num_inputs)] combined_inputs = imagen.deprecated.SeparatedComposite(min_separation=0,generators=inputs, bounds = retinal_bounds) if p.dataset in ["Nature", "VGR"]:

patch_orientation = 0.0 if p.dataset=="VGR" else numbergen.UniformRandom(lbound=-pi,ubound=pi,seed=p.input_seed+65)

if p.dataset=="Nature":
    image_filenames= [param.resolve_path("images/shouval/combined%02d.png"%(i+1)) for i in xrange(25)]
else:

    ordering = [1,25,8,24,22,23,4,15,17,13,5,20,18,14,3,6,12,10,21,7,9,2,16,11,19]
    image_filenames  = [param.resolve_path("images/VGR/combined%02d_VGR.png" % i) for i in ordering]

inputs=[pattern.image.FileImage(filename=f,
            scale=contrast_scale,
            size=10.0,
            x=numbergen.UniformRandom(lbound=-0.75,ubound=0.75,seed=p.input_seed+12),
            y=numbergen.UniformRandom(lbound=-0.75,ubound=0.75,seed=p.input_seed+36),
            orientation=patch_orientation)  for f in image_filenames]

combined_inputs = pattern.Selector(generators=inputs)

noise_ratio = (p.percent_noise / 100.0) disk_scale= 1.0 / (1.0 + noise_ratio) rand_scale= noise_ratio / (1.0 + noise_ratio)

disks_inputs=[topo.pattern.Composite(operator=numpy.add, scale=contrast_scale, generators=[topo.pattern.Disk( x=numbergen.UniformRandom(lbound=-2.125,ubound=2.125, seed=p.input_seed+12), y=numbergen.UniformRandom(lbound=-2.125,ubound=2.125, seed=p.input_seed+36), size=2.0, aspect_ratio=1.0, scale = disk_scale, bounds=sheet.BoundingBox(radius=1.125),smoothing=0.1), topo.pattern.random.UniformRandom(scale=rand_scale)])]

retina_inputs = topo.pattern.Selector(generators=disks_inputs)

if p.retinal_waves == 0: retina_inputs = combined_inputs else: topo.sim.schedule_command(p.retinal_waves, 'topo.sim["Retina"].set_input_generator(combined_inputs, push_existing=False)')

projection.CFProjection.cf_shape=pattern.Disk(smoothing=0.0) projection.CFProjection.response_fn=responsefn.optimized.CFPRF_DotProduct_opt() projection.CFProjection.learning_fn=learningfn.optimized.CFPLF_Hebbian_opt() projection.CFProjection.weights_output_fns=[transferfn.optimized.CFPOF_DivisiveNormalizeL1_opt()] projection.SharedWeightCFProjection.response_fn=responsefn.optimized.CFPRF_DotProduct_opt()

topo.sim['Retina']=sheet.GeneratorSheet(nominal_density=p.retina_density, input_generator=retina_inputs, period=1.0, phase=0.05, nominal_bounds=sheet.BoundingBox(points=((-(2.0_p.area/2.0+0.25+0.375+0.5),-(p.area/2.0+0.25+0.375+0.5)),((2.0_p.area/2.0+0.25+0.375+0.5),p.area/2.0+0.25+0.375+0.5))))

lgn_surroundg = pattern.Gaussian(size=0.25,aspect_ratio=1.0, output_fns=[transferfn.DivisiveNormalizeL1()])

for s in ['LGNOn','LGNOff']:

extra_kwargs = dict(tsettle=2,strict_tsettle=1) if p.gain_control else dict(tsettle=0,strict_tsettle=0)

topo.sim[s]=sheet.optimized.LISSOM_Opt(nominal_density=p.lgn_density,
               nominal_bounds=sheet.BoundingBox(points=((-(2.0*p.area/2.0+0.25+0.5),-(p.area/2.0+0.25+0.5)),((2.0*p.area/2.0+0.25+0.5),p.area/2.0+0.25+0.5))),
                output_fns=[transferfn.misc.HalfRectify()],
                measure_maps=False, **extra_kwargs)

if p.gain_control:

    topo.sim.connect(s,s,delay=0.05, name='LateralGC',
                dest_port=("Activity"),activity_group=(0.6, DivideWithConstant(c=0.11)),
                connection_type=projection.SharedWeightCFProjection,
                strength=0.6,weights_generator=lgn_surroundg,
                nominal_bounds_template=sheet.BoundingBox(radius=0.25))

learning_rate = 0.01 if p.homeostasis else 0.0

topo.sim["V1"] = sheet.lissom.LISSOM(nominal_density=p.cortex_density, tsettle=16, plastic=True, nominal_bounds=sheet.BoundingBox(points=((-2.0_p.area/2.0,-p.area/2.0),(2.0_p.area/2.0,p.area/2.0))), output_fns = [HomeostaticResponse(t_init=p.t_init, learning_rate=learning_rate)])

topo.sim["V1"].joint_norm_fn=topo.sheet.optimized.compute_joint_norm_totals_opt

centerg = pattern.Gaussian(size=0.07385, output_fns=[transferfn.DivisiveNormalizeL1()])

surroundg = pattern.Gaussian(size=0.29540, output_fns=[transferfn.DivisiveNormalizeL1()])

on_weights = pattern.Composite(generators=[centerg,surroundg],operator=numpy.subtract)

off_weights = pattern.Composite(generators=[surroundg,centerg],operator=numpy.subtract)

strength_factor = 6.0 topo.sim.connect( 'Retina','LGNOn',delay=0.05,strength=2.33*strength_factor, name='Afferent', connection_type=projection.SharedWeightCFProjection, nominal_bounds_template=sheet.BoundingBox(radius=0.375), weights_generator=on_weights)

topo.sim.connect( 'Retina','LGNOff',delay=0.05,strength=2.33*strength_factor, name='Afferent', connection_type=projection.SharedWeightCFProjection, nominal_bounds_template=sheet.BoundingBox(radius=0.375), weights_generator=off_weights)

"Center surround (difference-of-Gaussian) weights successfully generated"

LGN_V1_delay = 0.05 if p.gain_control else 0.10

topo.sim.connect( 'LGNOn','V1',delay=LGN_V1_delay,strength=p.aff_strength,name='LGNOnAfferent', dest_port=('Activity','JointNormalize','Afferent'), connection_type=projection.CFProjection,learning_rate=p.aff_lr, nominal_bounds_template=sheet.BoundingBox(radius=0.27083), weights_generator= pattern.random.GaussianCloud(gaussian_size=2*0.27083), learning_fn=learningfn.optimized.CFPLF_Hebbian_opt())

topo.sim.connect( 'LGNOff','V1',delay=LGN_V1_delay,strength=p.aff_strength,name='LGNOffAfferent', dest_port=('Activity','JointNormalize','Afferent'), connection_type=projection.CFProjection,learning_rate=p.aff_lr, nominal_bounds_template=sheet.BoundingBox(radius=0.27083), weights_generator= pattern.random.GaussianCloud(gaussian_size=2*0.27083), learning_fn=learningfn.optimized.CFPLF_Hebbian_opt())

"Afferent GaussianCloud weights successfully generated."

lateral_excitatory_weights = pattern.Gaussian(aspect_ratio=1.0, size=0.05)

topo.sim.connect( 'V1','V1',delay=0.05,strength=p.exc_strength,name='LateralExcitatory', connection_type=projection.CFProjection,learning_rate=p.exc_lr, nominal_bounds_template=sheet.BoundingBox(radius=0.104), weights_generator=lateral_excitatory_weights)

"Lateral excitatory weights successfully generated"

lateral_inhibitory_weights = pattern.random.GaussianCloud(gaussian_size=0.15)

topo.sim.connect( 'V1','V1',delay=0.05,strength=-1.0*p.inh_strength,name='LateralInhibitory', connection_type=projection.CFProjection,learning_rate=p.inh_lr, nominal_bounds_template=sheet.BoundingBox(radius=0.22917), weights_generator=lateral_inhibitory_weights)

"Lateral inhibitory weights successfully generated"

topo.sim.grid_layout([[None, 'V1', None], ['LGNOn', None, 'LGNOff'], [None, 'Retina', None]], xstart=150,item_scale=0.8)

import topo.analysis.featureresponses topo.analysis.featureresponses.FeatureMaps.selectivity_multiplier=2.0 topo.analysis.featureresponses.FeatureCurveCommand.contrasts= [100, 80, 60, 40, 20, 10]

from topo.command.analysis import measure_or_pref import save_data import param import os import errno

def make_sure_path_exists(path): try: os.makedirs(path) except OSError as exception: if exception.errno != errno.EEXIST: raise

path ="(not shown)" make_sure_path_exists(path) param.normalize_path.prefix = path

fname='s'+str(args.inputseed) filename = str(0).rjust(5, '0') for time in range(200): topo.sim.run(1000) measure_or_pref() p = topo.sim.V1.views.maps.OrientationPreference.last.data s = topo.sim.V1.views.maps.OrientationSelectivity.last.data n = time + 1 filename = str(n).rjust(5, '0') save_data.save(p,s,filename+'.txt',path+fname+'/')

windform2015 commented 7 years ago

Hi, could you help tell me that , are there any commands or functions in topographica to view the geometry of , V1, LGN and Retina ?

windform2015 commented 7 years ago

Now I find the model editor in GUI can do this. Still thanks.

jbednar commented 7 years ago

Right, the model editor us the current way to do that. Unfortunately, it is part of the Tk GUI that we no longer recommend for new projects, but if it works, great!