Source code

001package cnslab.cnsnetwork;
002import cnslab.cnsmath.*;
003import java.io.FileOutputStream;
004import java.io.ObjectOutputStream;
005import java.io.ObjectInputStream;
006import java.io.FileInputStream;
007import java.util.HashMap;
008import java.util.ArrayList;
009import java.util.Iterator;
010import java.util.Map;
011import java.util.Arrays;
012import java.util.TreeSet;
013/**
014 * A class of Utility functions.
015 * 
016 * @author Yi Dong
017 */
018public class FunUtil 
019{
020//      public static double minAccuracy=1e-8;
021        public static double minAccuracy=1e-10;
022
023        /**
024         * Binary search for host id, given an array of neuron index, which are the last neuron index for that nethost
025         * 
026         * @param endIndex
027         * @param num
028         * @return
029         */
030        public static int hostId(int [] endIndex, int num)
031        {
032                int begin=0;
033                int end=endIndex.length-1;
034                int mid;
035                if(num>endIndex[end])throw new RuntimeException("hostId, neuron index "+num+" is out of range"+endIndex[end]);
036
037                while(end-begin>1)
038                {
039                        mid=(begin+end)/2;
040                        if(endIndex[mid]==num)return mid;
041                        if(endIndex[mid]<num)
042                        {
043                                begin=mid;
044                        }
045                        if(endIndex[mid]>num)
046                        {
047                                end=mid;
048                        }
049                }
050                return (endIndex[begin]>=num ? begin:end );
051        }
052        /*
053
054        public static Neuron[] randomVSICLNetwork(int[] endIndex,int numSyn,int numSensory)
055        {
056                //int[] endIndex={ 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000};
057
058                int totalNum=endIndex[endIndex.length-1]+1;
059                //int numSyn=1000;
060                int numTasks=endIndex.length;
061                //int numSensory=10;
062
063                if(numSyn > totalNum) throw new RuntimeException("# of synapses are more than neurons");
064                Seed idum= new Seed(-3);
065
066                Neuron [] neurons;
067                Synapse [] synapses;
068                Branch [] branches; 
069
070                neurons = new Neuron[totalNum];
071
072                for(int neuIndex=0; neuIndex< totalNum; neuIndex++)
073                {
074                        HashMap<String, ArrayList<Synapse> > tmp_bran= new HashMap<String,ArrayList<Synapse> >();
075
076                        int number = (int)(Cnsran.ran2(idum)*numSyn);
077
078                        HashMap<Integer,Boolean> nore = new HashMap<Integer, Boolean>();
079                        HashMap<String,Double> delays = new HashMap<String,Double>();
080
081                        for(int i=0; i< number ;i++)
082                        {
083                                int x= (int)(Cnsran.ran2(idum)*3);
084                                double delay=0.002;
085                                if(x==0)delay=0.002;
086                                if(x==1)delay=0.003;
087                                if(x==2)delay=0.004;
088
089                                int targetId = (int)(Cnsran.ran2(idum)*totalNum);
090                                while(targetId==neuIndex || targetId < numSensory || nore.containsKey(targetId) )
091                                {
092                                        targetId = (int)(Cnsran.ran2(idum)*totalNum);
093                                }
094                                nore.put(targetId,true);
095
096                                int targetHostId=FunUtil.hostId(endIndex,targetId);
097
098                                if(tmp_bran.get(delay+";"+targetHostId)==null) 
099                                {
100                                        tmp_bran.put( delay+";"+targetHostId, (new ArrayList<Synapse>()));
101                                }
102                                delays.put(delay+";"+targetHostId,delay);
103
104                                if(neuIndex < numSensory)
105                                {
106                                        if(Cnsran.ran2(idum)<1)
107                                        {
108                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, 1e-10, 0)); //excitatory
109                                        }
110                                        else
111                                        {
112                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, -2e-13, 1));//inhibitory
113                                        }
114                                }
115                                else
116                                {
117                                        if(Cnsran.ran2(idum)<0.5)
118                                        {
119                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, 1e-10, 0)); //excitatory
120                                        }
121                                        else
122                                        {
123                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, -2e-11, 1));//inhibitory
124                                        }
125                                }
126                        }
127
128                        Iterator< Map.Entry<String, ArrayList<Synapse> >  > values = tmp_bran.entrySet().iterator();
129
130                        branches= new Branch[tmp_bran.size()];
131                        int iter=0;
132                        while(values.hasNext())
133                        {
134                                Map.Entry<String, ArrayList<Synapse> > entry= values.next();
135                                //System.out.println((values.next()).toArray()[0]);
136                                //System.out.println((values.next()).toArray().length);
137
138                                Synapse[] syns=(entry.getValue()).toArray(new Synapse [0]);
139                                //Arrays.sort(syns);
140                                branches[iter]= new Branch(syns, delays.get(entry.getKey())) ;
141                                iter++;
142                        }
143                        if(neuIndex<numSensory)
144                        {
145                                Seed idum2= new Seed(-neuIndex);
146                                neurons[neuIndex] = new BKPoissonNeuron(idum2,400.0,new Axon(branches));  // seed, background freq, axon 
147                        }
148                        else
149                        {
150                                neurons[neuIndex] = new VSICLIFNeuron(-0.07,0.01,-0.05,new double[] {0.0,0.0}, new Axon(branches));  // seed, background freq, axon 
151                        }
152                }
153                return neurons;
154        }
155        public static Neuron[] randomNetwork(int[] endIndex,int numSyn,int numSensory)
156        {
157                //int[] endIndex={ 1000,2000,3000,4000,5000,6000,7000,8000,9000,10000};
158
159                int totalNum=endIndex[endIndex.length-1]+1;
160                //int numSyn=1000;
161                int numTasks=endIndex.length;
162                //int numSensory=10;
163
164                if(numSyn > totalNum) throw new RuntimeException("# of synapses are more than neurons");
165                Seed idum= new Seed(-3);
166
167                Neuron [] neurons;
168                Synapse [] synapses;
169                Branch [] branches; 
170
171                neurons = new Neuron[totalNum];
172
173                for(int neuIndex=0; neuIndex< totalNum; neuIndex++)
174                {
175                        HashMap<String, ArrayList<Synapse> > tmp_bran= new HashMap<String,ArrayList<Synapse> >();
176
177                        int number = (int)(Cnsran.ran2(idum)*numSyn);
178
179                        HashMap<Integer,Boolean> nore = new HashMap<Integer, Boolean>();
180                        HashMap<String,Double> delays = new HashMap<String,Double>();
181
182                        for(int i=0; i< number ;i++)
183                        {
184                                int x= (int)(Cnsran.ran2(idum)*3);
185                                double delay=0.002;
186                                if(x==0)delay=0.002;
187                                if(x==1)delay=0.003;
188                                if(x==2)delay=0.004;
189
190                                int targetId = (int)(Cnsran.ran2(idum)*totalNum);
191                                while(targetId==neuIndex || targetId < numSensory || nore.containsKey(targetId) )
192                                {
193                                        targetId = (int)(Cnsran.ran2(idum)*totalNum);
194                                }
195                                nore.put(targetId,true);
196
197                                int targetHostId=FunUtil.hostId(endIndex,targetId);
198
199                                if(tmp_bran.get(delay+";"+targetHostId)==null) 
200                                {
201                                        tmp_bran.put( delay+";"+targetHostId, (new ArrayList<Synapse>()));
202                                }
203                                delays.put(delay+";"+targetHostId,delay);
204
205                                if(neuIndex < numSensory)
206                                {
207                                        if(Cnsran.ran2(idum)<1)
208                                        {
209                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, 1e-10, 0)); //excitatory
210                                        }
211                                        else
212                                        {
213                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, -2e-13, 1));//inhibitory
214                                        }
215                                }
216                                else
217                                {
218                                        if(Cnsran.ran2(idum)<0.5)
219                                        {
220                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, 1e-10, 0)); //excitatory
221                                        }
222                                        else
223                                        {
224                                                tmp_bran.get( delay+";"+targetHostId).add(new Synapse(targetId, -2e-11, 1));//inhibitory
225                                        }
226                                }
227
228                        }
229
230                        Iterator< Map.Entry<String, ArrayList<Synapse> >  > values = tmp_bran.entrySet().iterator();
231
232                        branches= new Branch[tmp_bran.size()];
233                        int iter=0;
234                        while(values.hasNext())
235                        {
236                                Map.Entry<String, ArrayList<Synapse> > entry= values.next();
237                                //System.out.println((values.next()).toArray()[0]);
238                                //System.out.println((values.next()).toArray().length);
239
240                                Synapse[] syns=(entry.getValue()).toArray(new Synapse [0]);
241                                //Arrays.sort(syns);
242                                branches[iter]= new Branch(syns, delays.get(entry.getKey())) ;
243                                iter++;
244                        }
245                        if(neuIndex<numSensory)
246                        {
247                                Seed idum2= new Seed(-neuIndex);
248                                neurons[neuIndex] = new BKPoissonNeuron(idum2,400.0,new Axon(branches));  // seed, background freq, axon 
249                        }
250                        else
251                        {
252                                neurons[neuIndex] = new SIFNeuron(-0.07,0.01,new Axon(branches));  // seed, background freq, axon 
253                        }
254                }
255                return neurons;
256        }
257        */
258
259        /**
260         *  left rotate matrix 90 degree 
261         * 
262         * @param matrix 
263         * @return mattrix array
264         */
265        public static double[][] lRotate90(double conn [][])
266        {
267                int m;
268                m=conn.length;
269                int n;
270                n=conn[0].length;
271                double [][] out= new double[n][m];
272
273                for(int i=0;i<n;i++)
274                {
275                        for(int j=0;j<m;j++)
276                        {
277                                out[i][j]=conn[j][n-i-1];
278                        }
279                }
280                return out;
281        }
282
283        /**
284         *  right rotate matrix 90 degree 
285         * 
286         * @param matrix 
287         * @return mattrix array
288         */
289        public static double[][] rRotate90(double conn [][])
290        {
291                int m;
292                m=conn.length;
293                int n;
294                n=conn[0].length;
295                double [][] out= new double[n][m];
296
297                for(int i=0;i<n;i++)
298                {
299                        for(int j=0;j<m;j++)
300                        {
301                                out[i][j]=conn[m-j-1][i];
302                        }
303                }
304                return out;
305        }
306
307        /**
308         * Compute the standard deviation of the vector
309         *  
310         * @param vector
311         * @return standard deviation
312         */
313        public static double sd(double []d)
314        {
315                return Math.sqrt(var(d));
316        }
317
318        /**
319         * Compute the standard deviation of the vector, starting at start
320         * @param starting index
321         * @param vector
322         * @return standard deviation
323         */
324        public static double sd(int start, double []d)
325        {
326                return Math.sqrt(var(start,d));
327        }
328
329        public static double sd(int start, int end, double []d)
330        {
331                return Math.sqrt(var(start,end,d));
332        }
333
334
335        /**
336         * compute the variance of the vector
337         * 
338         * @param  vector 
339         * @return variance
340         */
341        public static double var(double []d)
342        {
343                double tmp=0.0;
344                double mean=mean(d);
345                for( double a : d)
346                {
347                        tmp+=(a-mean)*(a-mean);
348                }
349                tmp/=(double)(d.length-1);
350                return tmp;     
351        }
352
353        /**
354         * compute the variance of the vector starting at index start
355         * 
356         * @param  vector 
357         * @return variance
358         */
359        public static double var(int start, double []d)
360        {
361                double tmp=0.0;
362                double mean=mean(start, d);
363                for(int i=start; i< d.length; i++)
364                {
365                        tmp+=(d[i]-mean)*(d[i]-mean);
366                }
367                tmp/=(double)(d.length-start-1);
368                return tmp;     
369        }
370
371        /**
372         * 
373         * compute the variance of the vector starting at index start, ending at end index inclusively
374         * @param start id
375         * @param end id
376         * @param vector
377         * @return variance
378         */
379        public static double var(int start, int end,  double []d)
380        {
381                double tmp=0.0;
382                double mean=mean(start, d);
383                for(int i=start; i <= end ; i++)
384                {
385                        tmp+=(d[i]-mean)*(d[i]-mean);
386                }
387                tmp/=(double)(end-start);
388                return tmp;     
389        }
390
391
392        /**
393         *  Same as var function, except mean is computed
394         * 
395         * @param d
396         * @return mean
397         */
398        public static double mean(double []d)
399        {
400                double tmp=0.0;
401                for( double a : d)
402                {
403                        tmp+=a;
404                }
405                tmp/=(double)d.length;
406                return tmp;
407        }
408
409        /**
410         *   Same as var function, except mean is computed
411         * 
412         * @param start id
413         * @param  mean
414         * @return
415         */
416        public static double mean(int start, double []d)
417        {
418                double tmp=0.0;
419                for( int i = start; i< d.length; i++)
420                {
421                        tmp+=d[i];
422                }
423                tmp/=(double)(d.length- start);
424                return tmp;
425        }
426
427        /**
428         * 
429         *   Same as var function, except mean is computed
430         * @param start
431         * @param end
432         * @param d
433         * @return
434         */
435        public static double mean(int start, int end, double []d)
436        {
437                double tmp=0.0;
438                for( int i = start; i <= end; i++)
439                {
440                        tmp+=d[i];
441                }
442                tmp/=(double)(end - start + 1);
443                return tmp;
444        }
445
446        /**
447         * compute the difference of vector one and vector two 
448         * 
449         * @param vector one
450         * @param vector two
451         * @return vector
452         */
453        public static double[] diff(double []d, double []e)
454        {
455                if(d.length!=e.length) throw new RuntimeException("array of same size is expected");
456                double [] tmp= new double [d.length];
457                for( int i = 0; i < d.length; i++)
458                {
459                        tmp[i] = d[i]-e[i];
460                }
461                return tmp;
462        }
463
464        @Deprecated
465        public static String Pout(double []coefficients)
466        {
467                String ans="";
468                if(coefficients.length>0)
469                        ans+=coefficients[0];
470
471                for(int a=1;a<coefficients.length;a++)
472                {
473                        if(coefficients[a]!=0.0)
474                                ans+=" + " + coefficients[a] + "*x^" + a;
475                }
476                return ans;
477        }
478
479        @Deprecated
480        public static String RPout(double []coefficients)
481        {
482                String ans="";
483                if(coefficients.length>0)
484                        ans+=coefficients[coefficients.length-1];
485
486                for(int a=coefficients.length-2; a>=0; a--)
487                {
488                        if(coefficients[a]!=0.0)
489                                ans+=" + " + coefficients[a] + "*x^" + (coefficients.length-a-1);
490                }
491                return ans;
492        }
493
494        @Deprecated
495        public static double evaluate(double []coefficients, double x)
496        {
497                double sum=0.0;
498                sum =  coefficients[coefficients.length-1];
499                for(int exponent=coefficients.length-2 ;exponent >= 0;exponent--) 
500                {
501                        sum = x*sum+coefficients[exponent];
502                }
503                return sum;
504        }
505
506        @Deprecated
507        public static double[] talorShiftOne(double []out) //g(x)=f(x+1)
508        {
509                //double [] out = new double[coefficients.length];
510                int last = out.length-1;
511//              out[last] = out[last];
512                for(int i=1; i< out.length; i++)
513                {
514                        double tmpcof = out[last-i];
515                        out[last-i] = out[last-i+1];
516
517                        double tmp= out[last];
518                        double tmp2= out[last-1];
519
520                        for(int j=1; j<i; j++)
521                        {
522                                out[last-j] = out[last-j]+tmp;
523                                tmp=tmp2;
524                                tmp2=out[last-j-1];
525                        }
526                        out[last] = out[last]+tmpcof;
527//                      for(int j=0;j<i+1;j++)
528//                              System.out.print(out[last-j]+" ");
529//                      System.out.println();
530                }
531                return out;
532        }
533
534        @Deprecated
535        public static double[] RtalorShiftOne(double []out) //g(x)=f(x+1)
536        {
537                //double [] out = new double[coefficients.length];
538                int last = out.length-1;
539//              out[last] = out[last];
540                for(int i=1; i< out.length; i++)
541                {
542                        double tmpcof = out[i];
543                        out[i] = out[i-1];
544
545                        double tmp= out[0];
546                        double tmp2= out[1];
547
548                        for(int j=1; j<i; j++)
549                        {
550                                out[j] = out[j]+tmp;
551                                tmp=tmp2;
552                                tmp2=out[j+1];
553                        }
554                        out[0] = out[0]+tmpcof;
555//                      for(int j=0;j<i+1;j++)
556//                              System.out.print(out[last-j]+" ");
557//                      System.out.println();
558                }
559                return out;
560        }
561
562        @Deprecated
563        public static double[] talorShiftNegOne(double []out) //g(x)=f(x+1)
564        {
565                //double [] out = new double[coefficients.length];
566                int last = out.length-1;
567//              out[last] = out[last];
568                for(int i=1; i< out.length; i++)
569                {
570                        double tmpcof = out[last-i];
571                        out[last-i] = out[last-i+1];
572
573                        double tmp= out[last];
574                        double tmp2= out[last-1];
575
576                        for(int j=1; j<i; j++)
577                        {
578                                out[last-j] = -out[last-j]+tmp;
579                                tmp=tmp2;
580                                tmp2=out[last-j-1];
581                        }
582                        out[last] = -out[last]+tmpcof;
583//                      for(int j=0;j<i+1;j++)
584//                              System.out.print(out[last-j]+" ");
585//                      System.out.println();
586                }
587                return out;
588        }
589
590        @Deprecated
591        public static double[] RtalorShiftNegOne(double []out) //g(x)=f(x+1)
592        {
593                //double [] out = new double[coefficients.length];
594                int last = out.length-1;
595//              out[last] = out[last];
596                for(int i=1; i< out.length; i++)
597                {
598                        double tmpcof = out[i];
599                        out[i] = out[i-1];
600
601                        double tmp= out[0];
602                        double tmp2= out[1];
603
604                        for(int j=1; j<i; j++)
605                        {
606                                out[j] = -out[j]+tmp;
607                                tmp=tmp2;
608                                tmp2=out[j+1];
609                        }
610                        out[0] = -out[0]+tmpcof;
611//                      for(int j=0;j<i+1;j++)
612//                              System.out.print(out[last-j]+" ");
613//                      System.out.println();
614                }
615                return out;
616        }
617
618        @Deprecated
619        public static double [] H(double []out, double n) //g(x)=f(nx)
620        {
621                if(n==1.0)return reverse(out);
622                int last = out.length-1;
623                for(int i=1; i< out.length; i++)
624                {
625                        double tmpcof = out[last-i];
626                        for(int j=i; j>= 1; j--)
627                        {
628                                out[last-j] = out[last-j+1]*n;
629                        }
630                        out[last] = tmpcof;
631                        /*
632                        for(int j=0;j<i+1;j++)
633                                System.out.print(out[j]+" ");
634                        System.out.println();
635                        */
636                }
637                return out;
638        }
639
640        @Deprecated
641        public static double [] RH(double []out, double n) //g(x)=f(nx)
642        {
643                if(n==1.0)return reverse(out);
644                int last = out.length-1;
645                for(int i=1; i< out.length; i++)
646                {
647                        double tmpcof = out[i];
648                        for(int j=i; j>= 1; j--)
649                        {
650                                out[j] = out[j-1]*n;
651                        }
652                        out[0] = tmpcof;
653                        /*
654                        for(int j=0;j<i+1;j++)
655                                System.out.print(out[j]+" ");
656                        System.out.println();
657                        */
658                }
659                return out;
660        }
661
662        @Deprecated
663        public static int count(double []coefficients) {
664                int count=0;
665                boolean lastNeg=false;
666                boolean lastPos=false;
667
668                if(coefficients[0]>0)
669                        lastPos=true;
670                else if(coefficients[0]<0)
671                        lastNeg=true;
672
673
674                for(int i=1; i< coefficients.length; i++)
675                {
676                        if(coefficients[i]>0&&lastNeg){
677                                count++;
678                                lastPos=true;
679                                lastNeg=false;
680                        }
681                        else if(coefficients[i]<0&&lastPos){
682                                count++;
683                                lastPos=false;
684                                lastNeg=true;
685                        }
686                        else if(coefficients[i]<0&&!lastPos){
687                                lastNeg=true;
688                        }
689                        else if(coefficients[i]>0&&!lastNeg){
690                                lastPos=true;
691                        }
692                }
693                return count;
694        }
695
696        @Deprecated
697        public static double [] scale(double []out, double n)
698        {
699                for(int i=0;i<out.length;i++)out[i]*=n;
700                return out;
701        }
702
703        @Deprecated
704        public static double pow2(int expo)
705        {
706                if(expo>0)
707                {
708                        return (double)(1<<expo);
709                }
710                else if(expo<0)
711                {
712                        return 1.0/((double)(1<<(-expo)));
713                }
714                else
715                {
716                        return 1.0;
717                }
718        }
719
720        @Deprecated
721        public static double[] df(double []coefficients)
722        {
723                double[] dPrimeCoefs=new double[coefficients.length-1];
724                for(int a=1;a<coefficients.length;a++)
725                {
726                                dPrimeCoefs[a-1]=coefficients[a]*(double)a;
727                }
728                return dPrimeCoefs;
729        }
730
731        @Deprecated
732        public static int desBound(double []coefficients) {
733                double [] out = new double[coefficients.length];
734                System.arraycopy(coefficients,0,out,0,coefficients.length);
735                return count(RtalorShiftOne(out));
736        }
737
738        @Deprecated
739        public static double smallestRoot(double []coefficients) { // return the smallerst Root in interval (0,1) or -1.0 for no root
740                if(count(coefficients)==0)return -1.0;
741//              System.out.println("root of "+Pout(coefficients));
742                TreeSet<KCelement> poolKC  = new TreeSet<KCelement>();
743                poolKC.add(new KCelement(0,0));
744                double [] Q = new double [coefficients.length];
745                System.arraycopy(coefficients,0,Q,0,coefficients.length);
746
747                KCelement first,old;
748                old = new KCelement(0,0);
749
750                while((first=poolKC.pollFirst())!=null)
751                {
752                        if(first.k==32){
753                                throw new RuntimeException("k:"+first.k+"c:"+first.c+" poly:"+Pout(coefficients));
754                                //continue; //it is a mutiply root
755                        }
756//                      System.out.println(first);
757//                      System.out.println(Pout(Q));
758                        if((int)(pow2(old.k-first.k)*(double)first.c-old.c)==1)
759                        {
760                                Q = (scale(RH(talorShiftOne(Q),pow2(old.k-first.k)),pow2(((coefficients.length-1)*(first.k-old.k)))));
761                        //      Q = Q.talorShiftOne().H(pow2(old.k-first.k)).scale(pow2(((coefficients.length-1)*(first.k-old.k))));
762                        }
763                        else if((int)(pow2(old.k-first.k)*(double)first.c-old.c)==0)
764                        {
765                                Q = reverse(scale(H(Q,pow2(old.k-first.k)),pow2(((coefficients.length-1)*(first.k-old.k)))));
766                        //      Q = Q.H(pow2(old.k-first.k)).scale(pow2(((coefficients.length-1)*(first.k-old.k))));
767                        }
768                        else
769                        {
770                                throw new RuntimeException("no shift ("+first.c+","+first.k+") old kc ("+old.c+","+old.k+")");
771                        }
772
773                        int countRoot = desBound(Q);
774//                      System.out.println("desBound count:"+countRoot);
775                        if(countRoot ==1)
776                        {
777//                              System.out.println("left:"+(double)first.c/(double)(1<<first.k)+" right:"+(double)(first.c+1)/(double)(1<<first.k));
778//
779                                return safeRootFind(coefficients, (double)first.c/(double)(1<<first.k) ,(double)(first.c+1)/(double)(1<<first.k));
780                        }
781                        else if (countRoot >1)
782                        {
783                                if(Math.abs(evaluate(Q,0.5))<minAccuracy)return (2*(double)first.c+1.0)/(double)(1<<first.k+1);
784                                poolKC.add(new KCelement(2*first.c,first.k+1));
785                                poolKC.add(new KCelement(2*first.c+1,first.k+1));
786//                              System.out.println("add two size "+poolKC.size());
787
788                        }
789                        old = first;
790                }
791//              System.out.println("not found");
792                return -1.0;
793        }
794
795        @Deprecated
796        public static double safeRootFind(double[] Q, double guess, double end)
797        {
798                double [] dQ= df(Q);
799                int MAXIT=100;
800                double f1,fh;
801                double xl,xh;
802                f1=evaluate(Q,guess);
803                fh=evaluate(Q,end);
804                if( (f1>0 && fh >0) || (f1<0 && fh <0) )throw new RuntimeException("not bracketed");
805                if(Math.abs(f1)<minAccuracy)return guess;
806                //if(Math.abs(fh)<minAccuracy)return end;
807                if(f1 <0.0)
808                {
809                        xl=guess;
810                        xh=end;
811                }
812                else
813                {
814                        xh=guess;
815                        xl=end;
816                }
817                double rts = 0.5*(xl+xh); 
818                double fvalue=evaluate(Q,rts);
819                double dfvalue=evaluate(dQ,rts);
820                double dxold = Math.abs(xh-xl);
821                double dx=dxold;
822                double temp;
823                for( int j=0; j< MAXIT; j++)
824                {
825                        if((((rts-xh)*dfvalue-fvalue)*((rts-xl)*dfvalue-fvalue)>0.0) || (Math.abs(2.0*fvalue) > Math.abs(dxold*dfvalue))) //out of range or slow convergent then bisection
826                        {
827                                dxold = dx;
828                                dx=0.5*(xh-xl);
829                                rts=xl+dx;
830                                if(xl == rts) return rts;
831                        }
832                        else //Newton method
833                        {
834                                dxold=dx;
835                                dx = fvalue/dfvalue;
836                                temp = rts;
837                                rts -= dx;
838                                if(temp == rts) return rts;
839                        }
840                        if( Math.abs(dx)<minAccuracy) return rts;
841                        fvalue=evaluate(Q,rts);
842                        dfvalue=evaluate(dQ,rts);
843                        if( fvalue<0)
844                        {
845                                xl=rts;
846                        }
847                        else
848                        {
849                                xh=rts;
850                        }
851                }
852                throw new RuntimeException("max iteration");
853        //      return 0.0;
854        }
855        @Deprecated
856        public static double [] reverse(double []out) {
857                int i = 0;
858                int j = out.length - 1;
859                double tmp;
860                while (j > i) {
861                        tmp = out[j];
862                        out[j] = out[i];
863                        out[i] = tmp;
864                        j--;
865                        i++;
866                }
867                return out;
868        }
869
870        @Deprecated
871        public static double largestRoot(double [] out) // return the largest Root in interval (0,1) or -1.0 for no root
872        {
873                if(count(out)==0)
874                {
875                        return -1.0;
876                }
877                else
878                {
879                        double smallest=smallestRoot(reverse(talorShiftNegOne(reflex(out))));
880                        if(smallest<0.0)
881                        {
882                                return -1.0;
883                        }
884                        else
885                        {
886                        return 1.0-  smallest;
887                        }
888                }
889        }
890
891        @Deprecated
892        public static double [] reflex(double [] out)
893        {
894                for( int i=1; i< out.length; i+=2)
895                {
896                        out[i]=-out[i];
897                }
898                return out;
899        }
900
901        @Deprecated
902        public static double quickNewton(double [] Q, double guess)
903        {
904                double [] dQ= df(Q);
905                double fvalue=evaluate(Q,guess);
906                double dfvalue=evaluate(dQ,guess);
907                if(dfvalue>0.0)return largestRoot(Q);
908//              System.out.println("f:"+fvalue);
909//              System.out.println("df:"+dfvalue);
910//              System.out.println("dx:"+(fvalue/dfvalue));
911                double dx;
912                do
913                {
914                dx = fvalue/dfvalue;
915                guess -= dx;
916                fvalue=evaluate(Q,guess);
917                dfvalue=evaluate(dQ,guess);
918//              if(dfvalue>0.0)return largestRoot(Q);
919                if(dfvalue>0.0)return -1.0;
920                /*
921//              for check accuracy
922                if(dfvalue>0.0)
923                {
924                        double result = largestRoot(Q);
925                        if(result != -1.0) throw new RuntimeException("exception is met for quickNetwon");
926                        return -1.0;
927                }
928                */
929                } while(Math.abs(fvalue)>minAccuracy);
930                return guess;
931        }
932        /**
933         * compute greatest common factor for  two numbers
934         * 
935         * @param a
936         * @return
937         */
938        public static int gcd(int a, int b)
939        {
940                int maxI,minI;
941                if(a>b)
942                {
943                        maxI = a;
944                        minI = b;
945                }
946                else
947                {
948                        maxI = b;
949                        minI = a;
950                }
951                int residule;
952                do
953                {
954                        residule = (maxI % minI);
955                        maxI = minI;
956                        minI = residule;
957
958                }while(residule != 0);
959                return maxI;
960        }
961
962        /**
963         * compute greatest common factor for an array of integers. 
964         * 
965         * @param a
966         * @return
967         */
968        public static int gcd(int []a)
969        {
970                if(a.length<2) throw new RuntimeException("array length no less than 2");
971                int cd = a[0];
972                for(int i=1; i<a.length; i++)
973                {
974                        cd = gcd(cd,a[i]);
975                }
976                return cd;
977        }
978
979        /**
980         * print array 
981         * 
982         * @param a
983         */
984        public static void printA (int []a)
985        {
986                for(int i=0;i<a.length;i++)
987                {
988                        System.out.print(a[i]+" ");
989                }
990        }
991
992        /**
993         * print array 
994         * 
995         * @param a
996         */
997        public static void printlnA (int []a)
998        {
999                for(int i=0;i<a.length;i++)
1000                {
1001                        System.out.print(a[i]+" ");
1002                }
1003                        System.out.println();
1004        }
1005
1006}