Source code

001    package cnslab.cnsnetwork;
002
003    import java.io.*;
004    import java.util.*;
005
006    import cnslab.image.imageIO;
007    import cnslab.image.Convolution;
008    import cnslab.cnsmath.*;
009    import edu.jhu.mb.ernst.model.ModelFactory;
010    import edu.jhu.mb.ernst.model.Synapse;
011
012    /***********************************************************************
013    * Stimulus class contains all types of stimulus.
014    * 
015    * The injected spikes are assumed to be inhomogeneous Poisson spikes,
016    * similar to BKPoissonNeuron.
017    * 
018    * @version
019    *   $Date: 2012-08-04 13:43:22 -0500 (Sat, 04 Aug 2012) $
020    *   $Rev: 104 $
021    *   $Author: croft $
022    * @author
023    *   Yi Dong
024    * @author
025    *   David Wallace Croft, M.Sc.
026    * @author
027    *   Jeremy Cohen
028    ***********************************************************************/
029    public final class  Stimulus
030      implements Serializable, Transmissible
031    ////////////////////////////////////////////////////////////////////////
032    ////////////////////////////////////////////////////////////////////////
033    {
034
035    private static final long
036      serialVersionUID = 0L;
037
038    // public variables
039
040    public Seed  seed;
041
042    public SimulatorParser  simulatorParser;
043
044    // private variables
045
046    private final ModelFactory  modelFactory;
047
048    private final double
049      strength,
050      onTime,
051      offTime,
052      freq,
053      decay; // the decay of the frequency
054
055    /** array of neuron gets input. */
056    private final List<Integer>
057      neuronIndexList;
058
059    private final List<Double>
060      relativeStrengthList;
061
062    // private ArrayList<Double> relativeStrH;
063
064    // private ArrayList<Double> relativeStrV;
065
066    private final int
067      type,
068      channel;
069
070    ////////////////////////////////////////////////////////////////////////
071    // constructor methods
072    ////////////////////////////////////////////////////////////////////////
073
074    public  Stimulus (
075      final ModelFactory     modelFactory,
076      final SimulatorParser  simulatorParser,
077      final double           strength,
078      final int              type,
079      final double           onTime,
080      final double           offTime,
081      final double           freq,
082      final double           decay,
083      final Seed             seed,
084      final int              channel )
085    ////////////////////////////////////////////////////////////////////////
086    {
087      this.modelFactory = modelFactory;
088
089      this.strength = strength;
090
091      this.type = type;
092
093      this.simulatorParser = simulatorParser;
094
095      this.onTime = onTime;
096
097      this.offTime = offTime;
098
099      this.freq = freq;
100
101      this.decay = decay;
102
103      this.seed = seed;
104
105      this.channel = channel;
106
107      neuronIndexList = new ArrayList<Integer> ( );
108
109      relativeStrengthList = new ArrayList<Double> ( );
110
111      // relativeStrH= new ArrayList<Double>();
112      // relativeStrV= new ArrayList<Double>();
113    }
114
115    ////////////////////////////////////////////////////////////////////////
116    // accessor methods
117    ////////////////////////////////////////////////////////////////////////
118
119    public int getChannel  ( ) { return channel; }
120
121    public double  getDecay    ( ) { return decay; }
122
123    public double  getFreq     ( ) { return freq; }
124
125    public double  getOffTime  ( ) { return offTime; }
126
127    public double  getOnTime   ( ) { return onTime; }
128
129    public double  getStrength ( ) { return strength; }    
130
131    public int     getType     ( ) { return type; }
132
133    ////////////////////////////////////////////////////////////////////////
134    ////////////////////////////////////////////////////////////////////////
135
136    public ArrayList<PInputEvent>  init ( final int  sourceID )
137    ////////////////////////////////////////////////////////////////////////
138    {  
139      final ArrayList<PInputEvent>  out = new ArrayList<PInputEvent> ( );
140
141      final int neuronIndexListSize = neuronIndexList.size ( );
142
143      for ( int  i = 0; i < neuronIndexListSize; i++ )
144      {
145        final int toNeuronIndex = neuronIndexList.get ( i ).intValue ( );
146
147        final double
148          relativeStrength = relativeStrengthList.get ( i ).doubleValue ( );
149
150        final double
151          synapseStrength = relativeStrength * strength; 
152
153        final Synapse  synapse = modelFactory.createSynapse (
154          toNeuronIndex,
155          ( byte ) type,
156          ( float ) synapseStrength );
157
158        final double  time = getNextTime ( onTime );
159
160        if ( time >= 0 )
161        {
162          out.add ( new PInputEvent ( time, synapse, sourceID ) );
163        }
164      }
165
166      return out;
167    }
168
169/*
170    public ArrayList<PInputEvent> initH(int sourceID)
171    ////////////////////////////////////////////////////////////////////////
172    {  
173    ArrayList<PInputEvent> out = new ArrayList<PInputEvent>();
174    Synapse syn;
175    for(int i = 0 ; i < neuIndex.size();i++)
176    {
177      syn = new Synapse(neuIndex.get(i),relativeStrH.get(i)*strength,type);
178      double time = getNextTime(onTime);
179      if(time>=0)
180      {
181      out.add(new PInputEvent(time,syn,sourceID));
182      }
183    }
184    return out;
185    }
186
187    public ArrayList<PInputEvent> initV(int sourceID)
188    ////////////////////////////////////////////////////////////////////////
189    {  
190    ArrayList<PInputEvent> out = new ArrayList<PInputEvent>();
191    Synapse syn;
192    for(int i = 0 ; i < neuIndex.size();i++)
193    {
194      syn = new Synapse(neuIndex.get(i),relativeStrV.get(i)*strength,type);
195      double time = getNextTime(onTime);
196      if(time>=0)
197      {
198      out.add(new PInputEvent(time,syn,sourceID));
199      }
200    }
201    return out;
202    }
203*/
204
205    /***********************************************************************
206    * Get the next time given current time.
207    * 
208    * Inhomogeneous Poisson process. 
209    ***********************************************************************/
210    public double  getNextTime ( final double  currTime )
211    ////////////////////////////////////////////////////////////////////////
212    {
213      double  time = currTime;
214
215      time += ( -Math.log ( Cnsran.ran2 ( seed ) ) / freq );
216
217      while ( Cnsran.ran2 ( seed )
218        > Math.exp ( -decay * ( time - onTime ) ) )
219      {
220        time += ( -Math.log ( Cnsran.ran2 ( seed ) ) / freq );
221
222        if ( time > offTime )
223        {
224          return -1.0;
225        }
226      }
227
228      if ( time > offTime )
229      {
230        return -1.0;
231      }
232
233      return time;
234    }
235
236    /***********************************************************************
237    * Stimulate a square of neurons.
238    * 
239    * @param x1
240    * @param y1
241    * @param x2
242    * @param y2
243    * @param pre
244    * @param suf
245    * @param hostId
246    ***********************************************************************/
247    public void  plotSquare (
248      final int     x1,
249      final int     y1,
250      final int     x2,
251      final int     y2,
252      final String  pre,
253      final String  suf,
254      final int     hostId )
255    ////////////////////////////////////////////////////////////////////////
256    {
257      int xTotal=0;
258
259      int yTotal=0;
260
261      int xstep=0;
262
263      int ystep=0;
264
265      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
266
267      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
268
269      if(xMul < 0 )
270      {           
271        xTotal=simulatorParser.xEdgeLength;
272        
273        xstep=-xMul;
274      }
275      else
276      {       
277        xTotal=simulatorParser.xEdgeLength*xMul;
278        
279        xstep=1;
280      }
281      
282      if(yMul < 0 )
283      {       
284        yTotal=simulatorParser.yEdgeLength;
285        
286        ystep=-yMul;
287      }
288      else
289      {       
290        yTotal=simulatorParser.yEdgeLength*yMul;
291        
292        ystep=1;
293      }
294
295      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
296      {
297        throw new RuntimeException("out of range error");
298      }
299    
300      // Draw four non-overlapping lines
301
302      final int xShift = x1 > x2 ? xstep : -xstep;
303
304      final int yShift = y1 > y2 ? ystep : -ystep;
305
306      plotLine( x1 , y1 , x1 , y2 + yShift, pre , suf , hostId );
307
308      plotLine( x1 , y2 , x2 + xShift, y2 , pre , suf , hostId );
309
310      plotLine( x2 , y2 , x2 , y1 - yShift, pre , suf , hostId );
311
312      plotLine( x2 , y1 , x1 - xShift, y1 , pre , suf , hostId );
313    }
314
315    /***********************************************************************
316    * Stimulate a filled square of neurons
317    ***********************************************************************/
318    public void  plotFilledSquare (
319      final int     x1,
320      final int     y1,
321      final int     x2,
322      final int     y2,
323      final String  pre,
324      final String  suf,
325      final int     hostId )
326    ////////////////////////////////////////////////////////////////////////
327    {
328      int xTotal=0;
329
330      int yTotal=0;
331
332      int xstep=0;
333
334      int ystep=0;
335
336      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
337
338      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
339
340      if(xMul < 0 )
341      {           
342        xTotal=simulatorParser.xEdgeLength;
343        xstep=-xMul;
344      }
345      else
346      {       
347        xTotal=simulatorParser.xEdgeLength*xMul;
348        xstep=1;
349      } 
350      if(yMul < 0 )
351      {       
352        yTotal=simulatorParser.yEdgeLength;
353        ystep=-yMul;
354      }
355      else
356      {       
357        yTotal=simulatorParser.yEdgeLength*yMul;
358        ystep=1;
359      }
360
361      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
362      {
363        throw new RuntimeException("out of range error");
364      }
365
366
367      for (int x = Math.min( x1 , x2 ); x <= Math.max( x1 , x2 ); x+= xstep)
368      {
369        plotLine( x , y1 , x, y2 , pre , suf , hostId );
370      }
371    }
372
373    /***********************************************************************
374    * Stimulate a line of neurons
375    ***********************************************************************/
376    public void  plotLine (
377      final int     x1,
378      final int     y1,
379      final int     x2,
380      final int     y2,
381      final String  pre,
382      final String  suf,
383      final int     hostId )
384    ////////////////////////////////////////////////////////////////////////
385    {
386      double length = Math.sqrt(
387        (double)(x2-x1)*(double)(x2-x1) +
388        (double)(y2-y1)*(double)(y2-y1));
389
390      if(length < 1e-5)  // point stimulus
391      {
392        int id;
393
394        if ( ( id = simulatorParser.layerStructure.cellmap_slow (
395          pre, suf, x1, y1 ) ) != -1 )
396        {
397          if(FunUtil.hostId(
398            simulatorParser.layerStructure.nodeEndIndices, id) ==  hostId)
399          {
400            relativeStrengthList.add(1.0);
401            neuronIndexList.add(id);
402          }
403        }
404        else
405        {
406          throw new RuntimeException("point is not in grid");
407        }
408        return;
409      }
410
411      int xTotal=0;
412
413      int yTotal=0;
414
415      int xstep=0;
416
417      int ystep=0;
418
419      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
420
421      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
422
423      if(xMul < 0 )
424      {           
425        xTotal=simulatorParser.xEdgeLength;
426        xstep=-xMul;
427      }
428      else
429      {       
430        xTotal=simulatorParser.xEdgeLength*xMul;
431        xstep=1;
432      } 
433      if(yMul < 0 )
434      {       
435        yTotal=simulatorParser.yEdgeLength;
436        ystep=-yMul;
437      }
438      else
439      {       
440        yTotal=simulatorParser.yEdgeLength*yMul;
441        ystep=1;
442      }
443
444      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
445      {
446        throw new RuntimeException("out of range error");
447      }
448
449      // Here we determine the neurons that are "hit" by the line.
450      // See the Wikipedia page on Bresenham's line algorithm.
451
452      // If steep, swap x's and y's.
453      boolean steep = Math.abs(y2 - y1) > Math.abs(x2 - x1);
454      int _xBegin, _xEnd, _yBegin, _yEnd;
455      if (!steep)
456      {
457        _xBegin = x1;
458        _xEnd = x2;
459        _yBegin = y1;
460        _yEnd = y2;
461      }
462      else
463      {
464        _xBegin = y1;
465        _xEnd = y2;
466        _yBegin = x1;
467        _yEnd = x2;
468      }
469
470      // If flipped, swap begins and ends.
471      boolean flipped = _xBegin > _xEnd;
472      int xBegin, xEnd, yBegin, yEnd;
473      if (!flipped)
474      {
475        xBegin = _xBegin;
476        xEnd = _xEnd;
477        yBegin = _yBegin;
478        yEnd = _yEnd;
479      }
480      else
481      {
482        xBegin = _xEnd;
483        xEnd = _xBegin;
484        yBegin = _yEnd;
485        yEnd = _yBegin;
486      }
487
488      if(xBegin % xstep !=0) xBegin+= xstep - xBegin % xstep ;
489      if(yBegin % ystep !=0) yBegin+= ystep - yBegin % ystep ;
490
491      double slope = (double) (yEnd - yBegin) / (double) (xEnd - xBegin);
492
493      for (int x = xBegin; x <= xEnd; x += xstep)
494      {  
495        int y =  (int) Math.round(yBegin + slope * (x - xBegin));
496        y = y - (y % ystep);
497
498        // If steep, swap x's and y's
499        int xNeuron = steep ? y : x;
500        int yNeuron = steep ? x : y;
501
502        int id = simulatorParser.layerStructure.cellmap_slow(
503          pre,suf,xNeuron, yNeuron);
504
505        if(id != -1)
506        {      
507          if(FunUtil.hostId(
508            simulatorParser.layerStructure.nodeEndIndices, id) ==  hostId)
509          {
510            // Calculate the shortest distance from the point
511            // (xNeuron, yNeuron) to the line between (x1, y1) and (x2, y2)
512            
513            double distance = Math.abs((x2 - x1) * (y1 - yNeuron)
514              - (x1 - xNeuron) * (y2 - y1)) / length;
515
516            if(channel<0) // If there is no specified channel
517            {
518              // The relative strength of the synapse between this stimulus
519              // and a neuron is a function [ distanceTune() ] 
520              // of the distance between the neuron and the line.
521              
522              relativeStrengthList.add( distanceTune ( distance ) );
523            }
524            else if(channel >=0 && channel <= 180)
525            {
526              double lineAngle = Math.acos((double)(x2-x1)/length);
527
528              double  channelAngle
529                = (double) channel * (Math.PI / (double) 180.0);
530
531              relativeStrengthList.add (
532                distanceTune (
533                  distance * Math.cos(lineAngle - channelAngle) ) );
534            }
535            else
536            {
537              throw new RuntimeException (
538                "channel range error (0-180 or <0)" );
539            }
540
541            neuronIndexList.add(id);
542          }
543        }
544      }
545    }
546
547    /***********************************************************************
548    *  Given image file, stimulate the layer of neurons accourding to image
549    *  intensity
550    ***********************************************************************/
551    public void  plotImg (
552      final String          pre,
553      final String          suf,
554      final String          filename,
555      final double [ ] [ ]  filter )
556    ////////////////////////////////////////////////////////////////////////
557    {
558      int [][] out;
559      imageIO imgIO = new imageIO();
560      out = imgIO.loadBWImage("imagelib/"+filename);
561      int [][] pic = Convolution.convolve_safe(out,filter);
562
563      int height = pic.length; 
564      int width = pic[0].length;
565
566      int xstep,ystep;
567      int xmulti,ymulti;
568      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
569      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
570
571      int xreso,yreso;
572
573      if(xmulti>0)
574      {       
575        xstep=1;
576        xreso=xmulti*simulatorParser.xEdgeLength;
577      }
578      else
579      {       
580        xstep=-xmulti;
581        xreso=simulatorParser.xEdgeLength;
582      }
583
584      if(ymulti>0)
585      {       
586        ystep=1;
587        yreso=ymulti*simulatorParser.yEdgeLength;
588      }
589      else
590      {       
591        ystep=-ymulti;
592        yreso=simulatorParser.yEdgeLength;
593      }
594
595      double heightRatio;
596      double widthRatio;
597
598      try {
599        heightRatio = (double)height/(double)(yreso/ystep);
600        widthRatio = (double)width/(double)(xreso/xstep);
601
602      }
603      catch(Exception ex) {
604        throw new RuntimeException("Pre:"+pre+" Suf:"+suf);
605      }
606
607
608      for( int x = 0  ; x <  xreso/xstep;  x++)
609      {
610        for( int y = 0  ; y < yreso/ystep;  y++)
611        {
612          // String neu = pre+","+x*xstep+","+ (yreso-y*ystep-1)+","+suf;
613
614          if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
615            [ ( int ) ( x * widthRatio  ) ] < 0 )
616          {
617            pic[(int)(y*heightRatio)+ystep-1][(int)(x*widthRatio)]=0;
618          }
619
620          if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
621            [ ( int ) ( x * widthRatio  ) ] > 255 )
622          {
623            pic[(int)(y*heightRatio)+ystep-1][(int)(x*widthRatio)]=255;
624          }
625
626          // double pixel = (double)(pic[(int)((y)*heightRatio
627          // +(int)(heightRatio/2.0))][(int)((x)*widthRatio
628          // +(int)(widthRatio/2.0))])/(double)256;
629
630          final double
631          pixel = ( double ) (
632            pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
633              [ ( int ) ( x * widthRatio  ) ] ) / ( double ) 256;
634
635          if ( pixel > 0.01 )
636          {
637            relativeStrengthList.add(pixel);
638
639            neuronIndexList.add (
640              simulatorParser.layerStructure.cellmap_slow (
641                pre,
642                suf,
643                x * xstep,
644                ( yreso - y * ystep - ystep ) ) );
645          }
646        }
647      }
648    }
649
650    /***********************************************************************
651    * Given image file, stimulate the layer of neurons accourding to reverse
652    * of image intensity.
653    ***********************************************************************/
654    public void  plotRImg (
655      final String  pre,
656      final String  suf,
657      final String  filename,
658      final double [ ] [ ]  filter ) //plot reverse img
659    ////////////////////////////////////////////////////////////////////////
660    {
661      //pas.xEdgeLength;
662      //pas.yEdgeLength;
663      int [][] out;
664      imageIO imgIO = new imageIO();
665      out = imgIO.loadBWImage("imagelib/"+filename);
666      for(int i=0;i<out.length;i++)
667      {
668        for(int j=0;j<out[0].length;j++)
669        {
670          out[i][j]=255-out[i][j];
671        }
672      }
673      int [][] pic = Convolution.convolve_safe(out,filter);
674      int height = pic.length; 
675      int width = pic[0].length;
676
677      int xstep,ystep;
678      int xmulti,ymulti;
679      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
680      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
681
682      int xreso,yreso;
683
684      if(xmulti>0)
685      {       
686        xstep=1;
687        xreso=xmulti*simulatorParser.xEdgeLength;
688      }
689      else
690      {       
691        xstep=-xmulti;
692        xreso=simulatorParser.xEdgeLength;
693      }
694
695      if(ymulti>0)
696      {       
697        ystep=1;
698        yreso=ymulti*simulatorParser.yEdgeLength;
699      }
700      else
701      {       
702        ystep=-ymulti;
703        yreso=simulatorParser.yEdgeLength;
704      }
705
706      double heightRatio = (double)height/(double)(yreso/ystep);
707      double widthRatio = (double)width/(double)(xreso/xstep);
708
709      for( int x = 0  ; x <  xreso/xstep;  x++)
710      {
711        for( int y = 0  ; y < yreso/ystep;  y++)
712        {
713          // String neu = pre+","+x+","+ (pas.yEdgeLength-y-1)+","+suf;
714
715          if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
716            [ ( int ) ( x * widthRatio  ) ] < 0 )
717          {
718            pic[(int)(y*heightRatio)+ystep-1][(int)(x*widthRatio)]=0;
719          }
720
721          if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
722            [ ( int ) ( x * widthRatio  ) ] > 255 )
723          {
724            pic[(int)(y*heightRatio)+ystep-1][(int)(x*widthRatio)]=255;
725          }
726
727          final double  pixel = ( double ) (
728            pic [ ( int ) ( y * heightRatio ) + ystep - 1 ]
729              [ ( int ) ( x * widthRatio  ) ] ) / ( double ) 256;
730
731          //System.out.println("Layer pre:"+pre+" suf:"+suf+" heightRaito"
732          // +heightRatio+"X:"+(int)(y*heightRatio+ystep-1)+"_"+(x*xstep)
733          // +"Y:"+(int)(x*widthRatio)+"_"+(yreso-y*ystep-ystep)+" rate:"
734          // +(double)(pic[(int)(y*heightRatio)+ystep-1]
735          // [(int)(x*widthRatio)])/(double)256+" fin X:"+x+" Y:"+y);
736
737          if ( pixel > 0.01 )
738          {
739            relativeStrengthList.add ( pixel );
740
741            neuronIndexList.add (
742              simulatorParser.layerStructure.cellmap_slow (
743                pre,
744                suf,
745                x*xstep,
746                ( yreso-y*ystep-ystep ) ) );
747          }
748        }
749      }
750    }
751
752    /***********************************************************************
753    * Stimulate the layer of neurons with a given matrix.
754    * 
755    * cenX, cenY defines the center position of the matrix in the layer
756    ***********************************************************************/
757    public void  plotMatrix (
758      final String          pre,
759      final String          suf,
760      final int             cenX,
761      final int             cenY,
762      final double [ ] [ ]  matrix,
763      final int             hostId )
764    ////////////////////////////////////////////////////////////////////////
765    {
766      final int [ ]  center = simulatorParser.layerStructure.closePoint (
767        cenX,
768        cenY,
769        pre,
770        suf );
771
772      int halfMaX = matrix.length/2;
773      int halfMaY = matrix[0].length/2;
774
775      int xstep,ystep;
776      int xmulti,ymulti;
777      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
778      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
779
780      int xreso,yreso;
781
782      if(xmulti>0)
783      {       
784        xstep=1;
785        xreso=xmulti*simulatorParser.xEdgeLength;
786      }
787      else
788      {       
789        xstep=-xmulti;
790        xreso=simulatorParser.xEdgeLength;
791      }
792
793      if(ymulti>0)
794      {       
795        ystep=1;
796        yreso=ymulti*simulatorParser.yEdgeLength;
797      }
798      else
799      {       
800        ystep=-ymulti;
801        yreso=simulatorParser.yEdgeLength;
802      }
803
804      for (
805        int  x = center [ 0 ] - halfMaX;
806        x < center [ 0 ]
807          + ( matrix.length % 2 == 0 ? halfMaX : halfMaX + 1 );
808        x++ )
809      {
810        for (
811          int  y = center [ 1 ] - halfMaY;
812          y < center [ 1 ]
813            + ( matrix [ 0 ].length % 2 == 0 ? halfMaY : halfMaY + 1 );
814          y++ )
815        {
816          int MaX=x-center[0]+halfMaX;
817          int MaY=y-center[1]+halfMaY;
818          int xx,yy;
819          xx=x%xreso;
820          yy=y%yreso;
821          if(xx<0)xx=xx+xreso;
822          if(yy<0)yy=yy+yreso;
823
824          final int  id = simulatorParser.layerStructure.cellmap_slow (
825            pre,
826            suf,
827            xx,
828            yy );
829
830          if(id != -1)
831          {
832            if ( FunUtil.hostId (
833              simulatorParser.layerStructure.nodeEndIndices,
834              id ) ==  hostId )
835            {
836              double pixel = matrix[MaX][MaY];
837
838              if(pixel!=0.0)
839              {
840                relativeStrengthList.add(pixel);
841                neuronIndexList.add(id);
842              }
843            }
844          }
845        }
846      }
847    }
848
849    public void  plotMatrix (
850      final String          pre,
851      final String          suf,
852      final int             cenX,
853      final int             cenY,
854      final double [ ] [ ]  matrix )
855    ////////////////////////////////////////////////////////////////////////
856    {
857      final int [ ]  center = simulatorParser.layerStructure.closePoint (
858        cenX,
859        cenY,
860        pre,
861        suf );
862
863      int halfMaX = matrix.length/2;
864
865      int halfMaY = matrix[0].length/2;
866
867      int xstep,ystep;
868      int xmulti,ymulti;
869      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
870      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
871
872      int xreso,yreso;
873
874      if(xmulti>0)
875      {       
876        xstep=1;
877        xreso=xmulti*simulatorParser.xEdgeLength;
878      }
879      else
880      {       
881        xstep=-xmulti;
882        xreso=simulatorParser.xEdgeLength;
883      }
884
885      if(ymulti>0)
886      {       
887        ystep=1;
888        yreso=ymulti*simulatorParser.yEdgeLength;
889      }
890      else
891      {       
892        ystep=-ymulti;
893        yreso=simulatorParser.yEdgeLength;
894      }
895
896      for (
897        int  x = center [ 0 ] - halfMaX;
898        x < center [ 0 ]
899          + ( matrix.length % 2 == 0 ? halfMaX : halfMaX + 1 );
900        x++ )
901      {
902        for (
903          int  y = center [ 1 ] - halfMaY;
904          y < center [ 1 ]
905            + ( matrix [ 0 ].length % 2 == 0 ? halfMaY : halfMaY + 1 );
906          y++ )
907        {
908          int MaX=x-center[0]+halfMaX;
909          int MaY=y-center[1]+halfMaY;
910          int xx,yy;
911          xx=x%xreso;
912          yy=y%yreso;
913          if(xx<0)xx=xx+xreso;
914          if(yy<0)yy=yy+yreso;
915
916          final int  id = simulatorParser.layerStructure.cellmap_slow (
917            pre,
918            suf,
919            xx,
920            yy );
921
922          if(id != -1)
923          {
924            double pixel = matrix[MaX][MaY];
925            if(pixel!=0.0)
926            {
927              relativeStrengthList.add(pixel);
928              neuronIndexList.add(id);
929            }
930          }
931        }
932      }
933    }
934
935    /***********************************************************************
936    * Stimulate the whole layer of neurons.
937    *
938    * Injecting noise to the whole layer.
939    ***********************************************************************/
940    public void  background (
941      final String  pre,
942      final String  suf )
943    ////////////////////////////////////////////////////////////////////////
944    {
945      int
946        xTotal = 0,
947        yTotal = 0,
948        xstep  = 0,
949        ystep  = 0,
950        xMul   = simulatorParser.layerStructure.getXmultiplier ( pre, suf ),
951        yMul   = simulatorParser.layerStructure.getYmultiplier ( pre, suf );
952
953      if ( xMul < 0 )
954      {           
955        xTotal = simulatorParser.xEdgeLength;
956
957        xstep = -xMul;
958      }
959      else
960      {       
961        xTotal = simulatorParser.xEdgeLength * xMul;
962
963        xstep = 1;
964      }
965
966      if ( yMul < 0 )
967      {       
968        yTotal = simulatorParser.yEdgeLength;
969
970        ystep = -yMul;
971      }
972      else
973      {       
974        yTotal = simulatorParser.yEdgeLength * yMul;
975
976        ystep = 1;
977      }
978
979      for ( int  x = 0; x < xTotal; x = x + xstep )
980      {
981        for ( int  y = 0; y < yTotal; y = y + ystep )
982        {
983          relativeStrengthList.add ( Double.valueOf ( 1.0 ) );
984
985          neuronIndexList.add (
986            simulatorParser.layerStructure.cellmap_slow (
987              pre,
988              suf,
989              x,
990              y ) );
991        }
992      }
993    }
994
995    /***********************************************************************
996    * Injecting noise to the whole layer.
997    ***********************************************************************/
998    public void  background (
999      final String  pre,
1000      final String  suf,
1001      final int     hostId )
1002    ////////////////////////////////////////////////////////////////////////
1003    {
1004      int xTotal=0;           
1005      int yTotal=0;
1006      int xstep=0;
1007      int ystep=0;
1008      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
1009      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
1010
1011
1012      if(xMul < 0 )
1013      {           
1014        xTotal=simulatorParser.xEdgeLength;
1015        xstep=-xMul;
1016      }
1017      else
1018      {       
1019        xTotal=simulatorParser.xEdgeLength*xMul;
1020        xstep=1;
1021      } 
1022      if(yMul < 0 )
1023      {       
1024        yTotal=simulatorParser.yEdgeLength;
1025        ystep=-yMul;
1026      }
1027      else
1028      {       
1029        yTotal=simulatorParser.yEdgeLength*yMul;
1030        ystep=1;
1031      }
1032
1033
1034      for( int x = 0  ; x < xTotal;  x=x+xstep)
1035      {
1036        for( int y = 0  ; y < yTotal ; y=y+ystep)
1037        {
1038          //        String neu = pre+","+x+","+ y+","+suf;
1039          int id = simulatorParser.layerStructure.cellmap_slow(pre,suf,x,y);
1040
1041          if ( FunUtil.hostId (
1042            simulatorParser.layerStructure.nodeEndIndices,
1043            id ) == hostId )
1044          {
1045            relativeStrengthList.add(1.0);
1046
1047            neuronIndexList.add (
1048              simulatorParser.layerStructure.cellmap_slow (
1049                pre,
1050                suf,
1051                x,
1052                y ) );
1053          }
1054        }
1055      }
1056    }
1057
1058    /***********************************************************************
1059    * Plot reverse image.
1060    ***********************************************************************/
1061    public void plotRImg (
1062      final String          pre,
1063      final String          suf,
1064      final String          filename,
1065      final double [ ] [ ]  filter,
1066      final double          gamma,
1067      final double          k,
1068      final int             hostId )
1069    ////////////////////////////////////////////////////////////////////////
1070    {
1071      //pas.xEdgeLength;
1072      //pas.yEdgeLength;
1073      int [][] out;
1074      imageIO imgIO = new imageIO();
1075
1076      out = imgIO.loadBWImage("imagelib/"+filename);
1077
1078      for(int i=0;i<out.length;i++)
1079      {
1080        for(int j=0;j<out[0].length;j++)
1081        {
1082          out[i][j]=255-out[i][j];
1083        }
1084      }
1085      int [][] pic = Convolution.convolve_safe(out,filter);
1086      int height = pic.length; 
1087      int width = pic[0].length;
1088
1089      int xstep,ystep;
1090      int xmulti,ymulti;
1091      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
1092      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
1093
1094      int xreso,yreso;
1095
1096      if(xmulti>0)
1097      {       
1098        xstep=1;
1099        xreso=xmulti*simulatorParser.xEdgeLength;
1100      }
1101      else
1102      {       
1103        xstep=-xmulti;
1104        xreso=simulatorParser.xEdgeLength;
1105      }
1106
1107      if(ymulti>0)
1108      {       
1109        ystep=1;
1110        yreso=ymulti*simulatorParser.yEdgeLength;
1111      }
1112      else
1113      {       
1114        ystep=-ymulti;
1115        yreso=simulatorParser.yEdgeLength;
1116      }
1117
1118      double heightRatio = (double)height/(double)(yreso/ystep);
1119      double widthRatio = (double)width/(double)(xreso/xstep);
1120
1121      for( int x = 0  ; x <  xreso/xstep;  x++)
1122      {
1123        for( int y = 0  ; y < yreso/ystep;  y++)
1124        {
1125          // String neu = pre+","+x+","+ (pas.yEdgeLength-y-1)+","+suf;
1126
1127          int id = simulatorParser.layerStructure.cellmap_slow (
1128            pre,
1129            suf,
1130            x*xstep,
1131            (yreso-y*ystep-ystep));
1132
1133          if ( FunUtil.hostId (
1134            simulatorParser.layerStructure.nodeEndIndices,
1135            id ) ==  hostId )
1136            // if(pas.ls.cellmap.containsKey(neu))
1137          {
1138            // if(pic[(int)(y*heightRatio+ystep-1)]
1139            //       [(int)(x*widthRatio)]<0)
1140            // pic[(int)(y*heightRatio+ystep-1)][(int)(x*widthRatio)]=0;
1141
1142            // if(pic[(int)(y*heightRatio+ystep-1)]
1143            //       [(int)(x*widthRatio)]>255)
1144            //   pic[(int)(y*heightRatio+ystep-1)][(int)(x*widthRatio)]=255;
1145
1146            // double pixel = (double)(
1147            // 255-pic[(int)(y*heightRatio)]
1148            // [(int)(x*widthRatio)])/(double)256;
1149
1150            double pixel = 0.0;
1151
1152            int ii,jj;
1153
1154            for(ii=0;ii<ystep;ii++)
1155            {
1156              for(jj=0;jj<xstep;jj++)
1157              {
1158                if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1159                  [ ( int ) ( x * widthRatio  ) + jj ] < 0 )
1160                {
1161                  pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1162                    [ ( int ) ( x * widthRatio  ) + jj ] = 0;
1163                }
1164
1165                if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1166                  [ ( int ) ( x * widthRatio  ) + jj ] > 255 )
1167                {
1168                  pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1169                    [ ( int ) ( x * widthRatio  ) + jj ] = 255;
1170                }
1171
1172                pixel += nonlinearFilter (
1173                  ( double ) ( pic
1174                    [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1175                      [ ( int ) ( x * widthRatio  ) + jj ] )
1176                      / ( double ) 256,
1177                      gamma,
1178                      k );
1179              }
1180            }
1181
1182            pixel = pixel /(double)(ystep*xstep);
1183
1184            // double pixel = (double)(pic[(int)(y*heightRatio)]
1185            // [(int)(x*widthRatio)])/(double)256;
1186
1187            if(pixel>0.01)
1188            {
1189              relativeStrengthList.add(pixel);
1190              neuronIndexList.add(id);
1191            }
1192          }
1193        }
1194      }
1195    }
1196
1197    /***********************************************************************
1198    * Non linear filter to simualte Retina 
1199    * 
1200    * @param x
1201    * @param gamma
1202    * @param k
1203    * @return
1204    ***********************************************************************/
1205    public double nonlinearFilter(double x, double gamma, double k)
1206    ////////////////////////////////////////////////////////////////////////
1207    {
1208      //    double gamma = 1.0;
1209      //    double k=0.02;
1210      if(gamma>0)
1211      {
1212        double var = Math.pow((x/k),gamma);
1213
1214        return var/(1+var);
1215      }
1216
1217      return x;
1218    }
1219
1220    public void  plotImg (
1221      final String          pre,
1222      final String          suf,
1223      final String          filename,
1224      final double [ ] [ ]  filter,
1225      final double          gamma,
1226      final double          k,
1227      final int             hostId )
1228    ////////////////////////////////////////////////////////////////////////
1229    {
1230      //pas.xEdgeLength;
1231      //pas.yEdgeLength;
1232      int [][] out;
1233      imageIO imgIO = new imageIO();
1234      out = imgIO.loadBWImage("imagelib/"+filename);
1235      int [][] pic = Convolution.convolve_safe(out,filter);
1236
1237      int height = pic.length; 
1238      int width = pic[0].length;
1239
1240      int xstep,ystep;
1241      int xmulti,ymulti;
1242      xmulti = simulatorParser.layerStructure.getXmultiplier(pre,suf);
1243      ymulti = simulatorParser.layerStructure.getYmultiplier(pre,suf);
1244
1245      int xreso,yreso;
1246
1247      if(xmulti>0)
1248      {       
1249        xstep=1;
1250        xreso=xmulti*simulatorParser.xEdgeLength;
1251      }
1252      else
1253      {       
1254        xstep=-xmulti;
1255        xreso=simulatorParser.xEdgeLength;
1256      }
1257
1258      if(ymulti>0)
1259      {       
1260        ystep=1;
1261        yreso=ymulti*simulatorParser.yEdgeLength;
1262      }
1263      else
1264      {       
1265        ystep=-ymulti;
1266        yreso=simulatorParser.yEdgeLength;
1267      }
1268
1269      double heightRatio = (double)height/(double)(yreso/ystep);
1270
1271      double widthRatio  = (double)width /(double)(xreso/xstep);
1272
1273      for( int x = 0  ; x <  xreso/xstep;  x++)
1274      {
1275        for( int y = 0  ; y < yreso/ystep;  y++)
1276        {
1277          // String neu = pre+","+x+","+ (pas.yEdgeLength-y-1)+","+suf;
1278
1279          int id = simulatorParser.layerStructure.cellmap_slow(
1280            pre,
1281            suf,
1282            x*xstep,
1283            (yreso-y*ystep-ystep));
1284
1285          if ( FunUtil.hostId (
1286            simulatorParser.layerStructure.nodeEndIndices,
1287            id ) ==  hostId )
1288            // if(pas.ls.cellmap.containsKey(neu))
1289          {
1290            // double pixel = (double)(255-pic[(int)(y*heightRatio)]
1291            // [(int)(x*widthRatio)])/(double)256;
1292
1293            // double pixel = nonlinearFilter((double)(
1294            // pic[(int)(y*heightRatio)+ystep-1]
1295            // [(int)(x*widthRatio)])/(double)256,gamma,k);
1296
1297            double pixel = 0.0;
1298
1299            int ii,jj;
1300
1301            for(ii=0;ii<ystep;ii++)
1302            {
1303              for(jj=0;jj<xstep;jj++)
1304              {
1305                if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1306                  [ ( int ) ( x * widthRatio  ) + jj ] < 0 )
1307                {
1308                  pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1309                    [ ( int ) ( x * widthRatio  ) + jj ] = 0;
1310                }
1311
1312                if ( pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1313                  [ ( int ) ( x * widthRatio  ) + jj ] > 255 )
1314                {
1315                  pic [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1316                    [ ( int ) ( x * widthRatio  ) + jj ] = 255;
1317                }
1318
1319                pixel += nonlinearFilter (
1320                  ( double ) ( pic
1321                    [ ( int ) ( y * heightRatio ) + ystep - 1 - ii ]
1322                      [ ( int ) ( x * widthRatio  ) +jj ] )
1323                      / ( double ) 256,
1324                      gamma,
1325                      k );
1326              }
1327            }
1328
1329            pixel = pixel / ( double ) ( ystep * xstep );
1330
1331            if ( pixel > 0.01 )
1332            {
1333              relativeStrengthList.add ( pixel );
1334
1335              neuronIndexList.add ( id );
1336            }
1337          }
1338        }
1339      }
1340    }
1341
1342    public void  plotSquare (
1343      final int     x1,
1344      final int     y1,
1345      final int     x2,
1346      final int     y2,
1347      final String  pre,
1348      final String  suf )
1349    ////////////////////////////////////////////////////////////////////////
1350    {
1351      int xTotal=0;
1352
1353      int yTotal=0;
1354
1355      int xstep=0;
1356
1357      int ystep=0;
1358
1359      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
1360
1361      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
1362
1363      if(xMul < 0 )
1364      {           
1365        xTotal=simulatorParser.xEdgeLength;
1366        xstep=-xMul;
1367      }
1368      else
1369      {       
1370        xTotal=simulatorParser.xEdgeLength*xMul;
1371        xstep=1;
1372      } 
1373      if(yMul < 0 )
1374      {       
1375        yTotal=simulatorParser.yEdgeLength;
1376        ystep=-yMul;
1377      }
1378      else
1379      {       
1380        yTotal=simulatorParser.yEdgeLength*yMul;
1381        ystep=1;
1382      }
1383
1384      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
1385      {
1386        throw new RuntimeException("out of range error");
1387      }
1388
1389      // Draw four non-overlapping lines
1390
1391      final int xShift = x1 > x2 ? xstep : -xstep;
1392
1393      final int yShift = y1 > y2 ? ystep : -ystep;
1394
1395      plotLine( x1 , y1 , x1 , y2 + yShift, pre , suf );
1396
1397      plotLine( x1 , y2 , x2 + xShift, y2 , pre , suf );
1398
1399      plotLine( x2 , y2 , x2 , y1 - yShift, pre , suf );
1400
1401      plotLine( x2 , y1 , x1 - xShift, y1 , pre , suf );
1402    }
1403
1404    public void  plotFilledSquare (
1405      final int     x1,
1406      final int     y1,
1407      final int     x2,
1408      final int     y2,
1409      final String  pre,
1410      final String  suf )
1411    ////////////////////////////////////////////////////////////////////////
1412    {
1413      int xTotal=0;
1414
1415      int yTotal=0;
1416
1417      int xstep=0;
1418
1419      int ystep=0;
1420
1421      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
1422
1423      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
1424
1425      if(xMul < 0 )
1426      {           
1427        xTotal=simulatorParser.xEdgeLength;
1428        xstep=-xMul;
1429      }
1430      else
1431      {       
1432        xTotal=simulatorParser.xEdgeLength*xMul;
1433        xstep=1;
1434      } 
1435      if(yMul < 0 )
1436      {       
1437        yTotal=simulatorParser.yEdgeLength;
1438        ystep=-yMul;
1439      }
1440      else
1441      {       
1442        yTotal=simulatorParser.yEdgeLength*yMul;
1443        ystep=1;
1444      }
1445
1446      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
1447      {
1448        throw new RuntimeException("out of range error");
1449      }
1450
1451
1452      for (int x = Math.min( x1 , x2 ); x <= Math.max( x1 , x2 ); x+= xstep)
1453      {
1454        plotLine( x , y1 , x, y2 , pre , suf );
1455      }
1456    }
1457
1458    public void  plotLine (
1459      final int     x1,
1460      final int     y1,
1461      final int     x2,
1462      final int     y2,
1463      final String  pre,
1464      final String  suf )
1465    ////////////////////////////////////////////////////////////////////////
1466    {    
1467      double length = Math.sqrt(
1468        (double)(x2-x1)*(double)(x2-x1)+
1469        (double)(y2-y1)*(double)(y2-y1));
1470
1471      if(length < 1e-5)  // point stimulus
1472      {
1473        int id;
1474        if ( ( id = simulatorParser.layerStructure.cellmap_slow (
1475          pre, suf, x1, y1 ) ) != -1 )
1476        {
1477          relativeStrengthList.add(1.0);
1478          neuronIndexList.add(id);
1479        }
1480        else
1481        {
1482          throw new RuntimeException("point is not in grid");
1483        }
1484        return;
1485      }
1486
1487      int xTotal=0;
1488
1489      int yTotal=0;
1490
1491      int xstep=0;
1492
1493      int ystep=0;
1494
1495      int xMul=simulatorParser.layerStructure.getXmultiplier(pre,suf);
1496
1497      int yMul=simulatorParser.layerStructure.getYmultiplier(pre,suf);
1498
1499      if(xMul < 0 )
1500      {           
1501        xTotal=simulatorParser.xEdgeLength;
1502        xstep=-xMul;
1503      }
1504      else
1505      {       
1506        xTotal=simulatorParser.xEdgeLength*xMul;
1507        xstep=1;
1508      } 
1509      if(yMul < 0 )
1510      {       
1511        yTotal=simulatorParser.yEdgeLength;
1512        ystep=-yMul;
1513      }
1514      else
1515      {       
1516        yTotal=simulatorParser.yEdgeLength*yMul;
1517        ystep=1;
1518      }
1519
1520      if ( x1> xTotal || x2 > xTotal || y1> yTotal || y2> yTotal )
1521      {
1522        throw new RuntimeException("out of range error");
1523      }
1524
1525      // Here we determine the neurons that are "hit" by the line.
1526      // See the Wikipedia page on Bresenham's line algorithm.
1527
1528      // If steep, swap x's and y's.
1529      boolean steep = Math.abs(y2 - y1) > Math.abs(x2 - x1);
1530      int _xBegin, _xEnd, _yBegin, _yEnd;
1531      if (!steep)
1532      {
1533        _xBegin = x1;
1534        _xEnd = x2;
1535        _yBegin = y1;
1536        _yEnd = y2;
1537      }
1538      else
1539      {
1540        _xBegin = y1;
1541        _xEnd = y2;
1542        _yBegin = x1;
1543        _yEnd = x2;
1544      }
1545
1546      // If flipped, swap begins and ends.
1547      boolean flipped = _xBegin > _xEnd;
1548      int xBegin, xEnd, yBegin, yEnd;
1549      if (!flipped)
1550      {
1551        xBegin = _xBegin;
1552        xEnd = _xEnd;
1553        yBegin = _yBegin;
1554        yEnd = _yEnd;
1555      }
1556      else
1557      {
1558        xBegin = _xEnd;
1559        xEnd = _xBegin;
1560        yBegin = _yEnd;
1561        yEnd = _yBegin;
1562      }
1563
1564      if(xBegin % xstep !=0) xBegin+= xstep - xBegin % xstep ;
1565      if(yBegin % ystep !=0) yBegin+= ystep - yBegin % ystep ;
1566
1567      double slope = (double) (yEnd - yBegin) / (double) (xEnd - xBegin);
1568
1569      for (int x = xBegin; x <= xEnd; x += xstep)
1570      {  
1571        int y =  (int) Math.round(yBegin + slope * (x - xBegin));
1572        y = y - (y % ystep);
1573
1574        // If steep, swap x's and y's
1575        int xNeuron = steep ? y : x;
1576        int yNeuron = steep ? x : y;
1577
1578        int id = simulatorParser.layerStructure.cellmap_slow(
1579          pre,suf,xNeuron, yNeuron);
1580
1581        if(id != -1)
1582        {        
1583          // Calculate the shortest distance from the
1584          // point (xNeuron, yNeuron) to the
1585          // line between (x1, y1) and (x2, y2)
1586          
1587          double distance = Math.abs(
1588            (x2 - x1) * (y1 - yNeuron) -
1589            (x1 - xNeuron) * (y2 - y1)) / length;
1590
1591          if(channel<0) // If there is no specified channel
1592          {
1593            // The relative strength of the synapse between this stimulus
1594            // and a neuron is a function [ distanceTune() ] 
1595            // of the distance between the neuron and the line.
1596            relativeStrengthList.add( distanceTune ( distance ) );
1597          }
1598          else if(channel >=0 && channel <= 180)
1599          {
1600            double lineAngle = Math.acos((double)(x2-x1)/length);
1601
1602            double
1603              channelAngle = (double) channel * (Math.PI / (double) 180.0);
1604
1605            relativeStrengthList.add (
1606              distanceTune (
1607                distance * Math.cos(lineAngle - channelAngle) ) );
1608          }
1609          else
1610          {
1611            throw new RuntimeException (
1612              "channel range error (0-180 or <0)" );
1613          }
1614
1615          neuronIndexList.add(id);
1616        }
1617      }
1618    }
1619
1620    private double distanceTune (double d)
1621    ////////////////////////////////////////////////////////////////////////
1622    {
1623      return 1 - Math.sqrt(2) * d;
1624    }
1625
1626    ////////////////////////////////////////////////////////////////////////
1627    ////////////////////////////////////////////////////////////////////////
1628    }