def test():
    import datasource
    import model
    import dataset

    f = Fitter()

    # define source of "experimental data"
    function = lambda s, i: (10.0*s/1.0)/(1.0 + s/1.0 + i/5.0)
    data_source = datasource.Generated_ScanDataSource(
        function,
        ['s', 'i'],
        'v',
        [scipy.logspace(-2, 2, 50), scipy.logspace(-2, 2, 10)],
        noise=0.3
    )

    # define model
    model = model.Equation_Model("Vmax*s/Ks/(1 + s/Ks + i/Ki)", ['s', 'i'])

    dataset.ScanDataSet('name', f, data_source, model)

    # specify the optimization algorithm, defaults to scipy_leastsq
    #~ alg = algorithm.scipy_leastsq()
    alg = algorithm.robust_biweight()
    f.setAlgorithm(alg)

    # specify the parameters to be fitted
    f.addParameter('Vmax', init=1.0, min=0, max=100)
    f.addParameter('Ks', init=1.0, min=0, max=10)
    f.addParameter('Ki', init=1.0, min=0, max=10)

    r = f.solve()
    r.writeOutput()

 def D2_setup_scan(self,min1,max1,step1,min2,max2,step2):
     self.P1min = min1
     self.P1max = max1
     self.P1steps = step1
     self.P2min = min2
     self.P2max = max2
     self.P2steps = step2
     self.P1range = scipy.logspace(scipy.log10(min1),scipy.log10(max1),step1)
     self.P2range = scipy.logspace(scipy.log10(min2),scipy.log10(max2),step2)

    def _getParamGridFixedEffectModel(self, G0, G1, link):
        if link == 'linear':
            param_grid = dict(alpha=0.5*sp.logspace(-5, 5, 20))
        elif link == 'logistic':
            param_grid = dict(C=sp.logspace(-5, 5, 20))
        else:
            assert False

        return param_grid

 def D3_setup_scan(self,min1,max1,step1,min2,max2,step2,min3,max3,step3):
     self.P1min = min1
     self.P1max = max1
     self.P1steps = step1
     self.P2min = min2
     self.P2max = max2
     self.P2steps = step2
     self.P3min = min3
     self.P3max = max3
     self.P3steps = step3
     self.P1range = scipy.logspace(scipy.log10(min1),scipy.log10(max1),step1)
     self.P2range = scipy.logspace(scipy.log10(min2),scipy.log10(max2),step2)
     self.P3range = scipy.logspace(scipy.log10(min3),scipy.log10(max3),step3)

def degree_distrib(net, deg_type="total", node_list=None, use_weights=True,
                   log=False, num_bins=30):
    '''
    Computing the degree distribution of a network.
    
    Parameters
    ----------
    net : :class:`~nngt.Graph` or subclass
        the network to analyze.
    deg_type : string, optional (default: "total")
        type of degree to consider ("in", "out", or "total").
    node_list : list or numpy.array of ints, optional (default: None)
        Restrict the distribution to a set of nodes (default: all nodes).
    use_weights : bool, optional (default: True)
        use weighted degrees (do not take the sign into account: all weights
        are positive).
    log : bool, optional (default: False)
        use log-spaced bins.
    
    Returns
    -------
    counts : :class:`numpy.array`
        number of nodes in each bin
    deg : :class:`numpy.array`
        bins
    '''
    ia_node_deg = net.get_degrees(node_list, deg_type, use_weights)
    ra_bins = sp.linspace(ia_node_deg.min(), ia_node_deg.max(), num_bins)
    if log:
        ra_bins = sp.logspace(sp.log10(sp.maximum(ia_node_deg.min(),1)),
                               sp.log10(ia_node_deg.max()), num_bins)
    counts,deg = sp.histogram(ia_node_deg, ra_bins)
    ia_indices = sp.argwhere(counts)
    return counts[ia_indices], deg[ia_indices]

def default_frequency_range(syslist):
    """Compute a reasonable default frequency range for frequency
    domain plots.

    Finds a reasonable default frequency range by examining the features
    (poles and zeros) of the systems in syslist.

    Parameters
    ----------
    syslist : list of Lti
        List of linear input/output systems (single system is OK)

    Returns
    -------
    omega : array
        Range of frequencies in rad/sec

    Examples
    --------
    >>> from matlab import ss
    >>> sys = ss("1. -2; 3. -4", "5.; 7", "6. 8", "9.")
    >>> omega = default_frequency_range(sys)
    """
    # This code looks at the poles and zeros of all of the systems that
    # we are plotting and sets the frequency range to be one decade above
    # and below the min and max feature frequencies, rounded to the nearest
    # integer.  It excludes poles and zeros at the origin.  If no features
    # are found, it turns logspace(-1, 1)
    
    # Find the list of all poles and zeros in the systems
    features = np.array(())
    
    # detect if single sys passed by checking if it is sequence-like
    if (not getattr(syslist, '__iter__', False)):
        syslist = (syslist,)

    for sys in syslist:
        try:
            # Add new features to the list
            features = np.concatenate((features, np.abs(sys.pole())))
            features = np.concatenate((features, np.abs(sys.zero())))
        except:
            pass

    # Get rid of poles and zeros at the origin
    features = features[features != 0];

    # Make sure there is at least one point in the range
    if (features.shape[0] == 0): features = [1];

    # Take the log of the features
    features = np.log10(features)

    #! TODO: Add a check in discrete case to make sure we don't get aliasing
                        
    # Set the range to be an order of magnitude beyond any features
    omega = sp.logspace(np.floor(np.min(features))-1, 
                        np.ceil(np.max(features))+1)   

    return omega

 def fit(self, kk=None):
     """
     Fit Fourier spectrum with the function set at class instantination
     ==> NB: fitting is done in logarithmic coordinates
     and fills plotting arrays with data
     --------
     Options:
     --------
     kk
        (k1,k2) <None> spectral interval for function fitting
        by default interval [ kk[1], kk[imax__kk] ] will be fitted
        ==> i.e. k=0 is excluded
     """
     # fitting interval
     if kk:
         ik_min=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[0]).nonzero()[0][-1]
         ik_max=(self.fft_data.kk[1:self.fft_data.imax__kk]<=kk[1]).nonzero()[0][-1]
     else:
         ik_min=1;
         ik_max=self.fft_data.imax__kk
     # do fitting
     self.__popt,self.__pcov = scipy.optimize.curve_fit(self.__func_fit,
                                                        scipy.log(self.fft_data.kk[ik_min:ik_max]),
                                                        scipy.log(self.fft_data.Ik[ik_min:ik_max]) )
     # boundaries of fitted interval
     self.kmin = self.fft_data.kk[ik_min]
     self.kmax = self.fft_data.kk[ik_max]
     # fill plot arrays <===============
     self.kk_plot=scipy.logspace( scipy.log10(self.kmin),
                                  scipy.log10(self.kmax),
                                  self.nk_plot )
     self.Ik_plot=self.fitting_function(self.kk_plot)

def compute_rls(data_file, output_filename, rls_type=DEFAULT_RLS, save_out=DEFAULT_SAVE):
    """
    data file contains sampels and labels saved as a mat file
    with respective keywords. Make your own data wrapper if needed
    """

    if path.splitext(output_filename)[-1] != ".mat":
        output_filename += ".mat"

    if path.splitext(data_file)[-1] != ".mat":
        raise ValueError, "mat file needed"

    data = loadmat(data_file)
    X = data["samples"]
    Y = data["labels"]
    lambdas = sp.logspace(-6, 6, 30)

    if rls_type.lower() == "linear":
        w, loos = lrlsloo(X, Y, lambdas)
    elif rls_type.lower() == "nonlinear":
        w, loos = rlsloo(X, Y, lambdas)
    else:
        print "ERROR: specify linear or nonlinear"

    if save_out:
        out_data = {"weights": w, "loos": loos}
        savemat(out_fname, out_data)

    def run(self, npts=25, inv_points=None, access_limited=True, **kwargs):
        r"""
        Parameters
        ----------
        npts : int (default = 25)
            The number of pressure points to apply.  The list of pressures
            is logarithmically spaced between the lowest and highest throat
            entry pressures in the network.

        inv_points : array_like, optional
            A list of specific pressure point(s) to apply.

        """
        if 'inlets' in kwargs.keys():
            logger.info('Inlets recieved, passing to set_inlets')
            self.set_inlets(pores=kwargs['inlets'])
        if 'outlets' in kwargs.keys():
            logger.info('Outlets recieved, passing to set_outlets')
            self.set_outlets(pores=kwargs['outlets'])
        self._AL = access_limited
        if inv_points is None:
            logger.info('Generating list of invasion pressures')
            min_p = sp.amin(self['throat.entry_pressure']) * 0.98  # nudge down
            max_p = sp.amax(self['throat.entry_pressure']) * 1.02  # bump up
            inv_points = sp.logspace(sp.log10(min_p),
                                     sp.log10(max_p),
                                     npts)

        self._npts = sp.size(inv_points)
        # Execute calculation
        self._do_outer_iteration_stage(inv_points)

def create_grid(r_in, r_out, nshell, space = 'powerlaw1', end = True):
    # function to create grid
    if space == 'log10':
        from scipy import log10, logspace
        # get the exponent of the start- and
        # stop-radius in input units
        start = [log10(r_in), 0][r_in == 0]
        stop = log10(r_out)
        radii = logspace(start, stop, num=nshell, endpoint=end)
    elif space == "powerlaw1":
        from scipy import arange
        radii = r_in * (r_out/r_in)**(arange(nshell)/(nshell - 1.0))
    elif space == 'linear':
        from scipy import linspace
        # linearly spaced grid
        radii = linspace(r_in, r_out, num=nshell, endpoint=end)
    elif space == 'powerlaw2':
        from scipy import linspace
        # first check if coefficients to the power-law was given
        #~ if 'exp' in kwargs:
            #~ p_exp = kwargs['exp']
        #~ else: # if not, set it to 2, i.e. r^2
            #~ p_exp = 2
        radii = r_in + (r_out - r_in)*(linspace(r_in, r_out, num=nshell, endpoint=end)/(r_out))**2
        #pr_int('Not implemented yet.')
        #raise ParError(spaced)
    else:
        raise Exception(space)
    return radii

	def add_axis(self, param, start, stop, steps, logspace=False):
		"""Add a parameter discretization axis to the the grid

		Arguments
		---------
		param : string
			The name of the model parameter.
		start : float
			The starting value of the model parameter.
		stop : float
			The ending value of the model parameter.
		steps : integer
			The number of steps to insert between the start and stop.
		logspace : boolean
			Space the steps logarithmically?

		Returns
		-------
		None
		"""

		assert param in self.sim.get_model_params().keys()
		
		if logspace:
			self._grid_pts.append( scipy.logspace(start,stop,steps) )
		else:
			self._grid_pts.append( scipy.linspace(start,stop,steps) )
		self._grid_size *= len(self._grid_pts[-1])
		self._grid_order.append(param)

 def _make_forces(self, extreme_forces, num):
     a = extreme_forces[0]
     b = extreme_forces[1]
     c = extreme_forces[2]
     #print a, b, c, num
     forces = numpy.append(
         b -( scipy.logspace(scipy.log10(a), scipy.log10(b),
         num=num/2+1) - a),
         scipy.logspace(scipy.log10(b), scipy.log10(c),
         num=num-num/2)[1:]
         )
     forces = numpy.array( sorted(list(set(forces))) )
     if len(forces) != num:
         Exception("problem with forces: length ", len(forces))
     #print forces
     #for i in range(len(forces)):
     #    f = forces[i]
     #    rand = random.uniform( 0.99*f, 1.01*f)
     #    forces[i] = rand
     return forces

    def plot(self, indice):
        self.diagrama.ax.clear()
        self.diagrama.ax.set_xlim(0.1, 10)
        self.diagrama.ax.set_ylim(0, 1)
        self.diagrama.ax.set_xscale("log")

        self.diagrama.ax.set_title(QtWidgets.QApplication.translate(
            "pychemqt", "Heat Transfer Temperature Effectiveness"), size='12')
        self.diagrama.ax.set_xlabel("NTU", size='12')
        self.diagrama.ax.set_ylabel("P", size='14')
        self.diagrama.ax.set_xticklabels(["0.1", "1.0", "10"])
        xticklabels = ["0.2", "0.3", "", "0.5", "", "0.7", "", "", "2.0",
                       "3.0", "", "5.0", "", "7.0", "", ""]
        self.diagrama.ax.set_xticklabels(xticklabels, minor=True)

        flujo = self.flujo[indice][1]
        self.mixed.setVisible(flujo == "CrFSMix")
        kwargs = {}
        if flujo == "CrFSMix":
            kwargs["mixed"] = str(self.mixed.currentText())

        R = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.2, 1.4, 1.6, 1.8,
             2., 2.5, 3., 4., 6., 8., 10., 15.]

        NTU = logspace(-1.5, 1, 100)
        for ri in R:
            e = [0]+[TemperatureEffectiveness(N, ri, flujo, **kwargs) for N in NTU[1:]]
            self.diagrama.plot(NTU, e, "k")
            self.diagrama.ax.annotate(" R=%0.1f" % ri, (NTU[-1], e[-1]),
                                      size="medium", ha="left", va="center")

#        F=[0.3]
#        for f in F:
#            p=[]
#            NTU=[]
#            for r in R:
#                func=lambda P: CorrectionFactor(P, r, flujo)-f
#                print func(0.)
#                pi=fsolve(func, 0.2)
#                p.append(pi)
#                NTU.append(NTU_fPR(p, r, flujo))

#            p=[fsolve(lambda P: CorrectionFactor(P, r, flujo)-f, 0.5)[0] for r in R]
#            NTU=[NTU_fPR(pi, ri) for pi, ni in zip(p, R)]
#            self.diagrama.plot(NTU, p, "--")

        self.diagrama.draw()

        if flujo == "CrFSMix" and self.mixed.currentIndex():
            img = image.imread('images/equation/%s2.png' % flujo)
        else:
            img = image.imread('images/equation/%s.png' % flujo)
        self.image.set_data(img)
        self.refixImage()

    def run(self, npts=25, inv_pressures=None):
        r"""
        Run the algorithm for specified number of points or at given capillary
        pressures.

        Parameters
        ----------
        npts : scalar
            The number of points to obtain on the curve.  The points are
            automatically selected to span the range of capillary pressures
            using a logarithmic spacing (more points are lower capillary
            pressure values).

        inv_pressures : array_like
            A list of capillary pressures to apply. List should contain
            increasing and unique values.
        """
        # If no invasion points are given then generate some
        if inv_pressures is None:
            logger.info('Generating list of invasion pressures')
            min_p = sp.amin(self['throat.entry_pressure']) * 0.98  # nudge down
            max_p = sp.amax(self['throat.entry_pressure']) * 1.02  # bump up
            inv_points = sp.logspace(sp.log10(min_p),
                                     sp.log10(max_p),
                                     npts)
        else:
            # Make sure the given invastion points are sensible
            inv_points = sp.unique(inv_pressures)
        self._inv_points = inv_points

        # Ensure inlets are set
        if sp.sum(self['pore.inlets']) == 0:
            raise Exception('Inlet pores have not been specified')

        # Ensure outlet pores are set if trapping is enabled
        if self._trapping:
            if sp.sum(self['pore.outlets']) == 0:
                raise Exception('Outlet pores have not been specified')

        # Generate curve from points
        for inv_val in self._inv_points:
            # Apply one applied pressure and determine invaded pores
            logger.info('Applying capillary pressure: ' + str(inv_val))
            self._apply_percolation(inv_val)
            if self._trapping:
                logger.info('Checking for trapping')
                self._check_trapping(inv_val)

        # Find invasion sequence values (to correspond with IP algorithm)
        Pinv = self['pore.inv_Pc']
        self['pore.inv_seq'] = sp.searchsorted(sp.unique(Pinv), Pinv)
        Tinv = self['throat.inv_Pc']
        self['throat.inv_seq'] = sp.searchsorted(sp.unique(Tinv), Tinv)

 def setScanJobs(self, start, end, intervals, job, log=False):
     """Splits a range into a number of jobs with intervals"""
     assert intervals >= 1, '\n* Minimum of 1 interval'
     if log:
         kpoints = scipy.logspace(scipy.log10(start), scipy.log10(end), intervals+1)
     else:
         kpoints = scipy.linspace(start, end, intervals+1)
     self.job_list = []
     for p in range(len(kpoints)-1):
         job2 = job % (kpoints[p], kpoints[p+1])
         self.job_list.append(job2)
         print(job2)

    def runContinuation(self,parameter,low,high,density,par3d=None,logrange=True,runQuiet=True):
        """
        Run the continuation using the following parameters:

        Args:

        - parameter = str(the parameter to be scanned)
        - low = float(lower bound)
        - high = float(upper bound)
        - density = int(the number of initial points)
        - par3d = float(extra 3d parameter to insert into the output array) this parameter is not set ONLY used in output
        - logrange = boolean [default = True], if True generate the result using logspace(log10(low), log10(high), density) otherwise use a linear range
        - runQuiet = boolean [default = True], if True do not display intermediate results to screen, disable for debugging

        After running the continuation the results are stored in numpy arrays

        - mod.res_idx  = scan parameter values (and optionally par3d)
        - mod.res_metab = steady-state species concentrations
        - mod.res_flux = steady-state flux values

        """

        self.pitcon_scan_density = density
        self.pitcon_scan_parameter = parameter
        self.pitcon_scan_parameter_3d = par3d
        if logrange:
            self.pitcon_range_low = scipy.log10(low)
            self.pitcon_range_high = scipy.log10(high)
            self.model.pitcon_par_space = scipy.logspace(self.pitcon_range_low,self.pitcon_range_high,self.pitcon_scan_density)
        else:
            self.pitcon_range_low = low
            self.pitcon_range_high = high
            self.model.pitcon_par_space = scipy.linspace(self.pitcon_range_low,self.pitcon_range_high,self.pitcon_scan_density)

        self.model.pitcon_flux_gen = 1
        if runQuiet:
            self.model.SetQuiet()
        else:
            self.model.SetLoud()

        if self.pitcon_scan_parameter_3d != None:
            self.pitcon_res = self.model.PITCON(self.pitcon_scan_parameter, self.pitcon_scan_parameter_3d)
            self.res_idx = self.pitcon_res[:,:2]
            self.res_metab = self.pitcon_res[:,2:len(self.model.species)+2:]
            self.res_flux = self.pitcon_res[:,len(self.model.species)+2:]
        else:
            self.pitcon_res = self.model.PITCON(self.pitcon_scan_parameter)
            self.res_idx = self.pitcon_res[:,0]
            self.res_idx = self.res_idx.reshape(self.res_idx.shape[0],1)
            self.res_metab = self.pitcon_res[:,1:len(self.model.species)+1]
            self.res_flux = self.pitcon_res[:,len(self.model.species)+1:]

        print '\n\tContinuation complete\n'

def run() :

    # parameter dictionary
    p = {}                 
    p['number_batches'] = 3       
    p['leakage_penalty'] = 0.04
    p['assembly_width'] = 21.5036
    p['assembly_power'] = 3.4/193
    p['active_height'] = 366.0
    p['fuel_radius'] = 0.4096
    p['cladding_inner_radius'] = 0.4180
    p['cladding_outer_radius'] = 0.4750
    p['number_pins'] = 264
    p['power_share'] = 'reactivity'
        
    T_F = 900*sp.ones(p['number_batches'])    # batch fuel temperatures (K)
    T_C = 580*sp.ones(p['number_batches'])    # batch moderator temperatures (K)
    
    num_thick = 20
    #thick = sp.linspace(0.0, 500, num_thick)
    thick = sp.logspace(-1, sp.log10(5*10**2), num_thick)
    T_F_FeCrAl = sp.zeros((3, num_thick))
    T_C_FeCrAl = sp.zeros((3, num_thick))
    T_F_SiC = sp.zeros((3, num_thick))
    T_C_SiC = sp.zeros((3, num_thick))
    PPF_FeCrAl = sp.zeros((3, num_thick))
    PPF_SiC = sp.zeros((3, num_thick))
    
    solver = NRM(p, rho=rho, m2=m2, k_cladding=k_cladding)

    
    for i in range(num_thick) :    
        p['t_fecral'] = thick[i]
        p['t_sic'] = 0.0
        B, ppf, T_F, T_C = solver.solve(T_F, T_C)
        T_F_FeCrAl[:, i] = T_F[:]
        T_C_FeCrAl[:, i] = T_C[:]
        PPF_FeCrAl[:, i] = ppf[:]
    
        p['t_fecral'] = 0.0
        p['t_sic'] = thick[i]
        B, ppf, T_F, T_C = solver.solve(T_F, T_C)
        T_F_SiC[:, i] = T_F[:]
        T_C_SiC[:, i] = T_C[:]
        PPF_SiC[:, i] = ppf[:]
    
    pickle.dump({'thick': thick,
                 'T_F_FeCrAl': T_F_FeCrAl,
                 'T_C_FeCrAl': T_C_FeCrAl,
                 'PPF_FeCrAl': PPF_FeCrAl,
                 'T_F_SiC': T_F_SiC,
                 'T_C_SiC': T_C_SiC,
                 'PPF_SiC': PPF_SiC}, open('example_3.p', 'w'))

 def define_bins(self, **kwargs):
     r"""
     This defines the bins for a logscaled histogram
     """
     self.data_vector.sort()
     sf = self.args['scale_fact']
     num_bins = int(sp.logn(sf, self.data_vector[-1]) + 1)
     #
     # generating initial bins from 1 - sf**num_bins
     low = list(sp.logspace(0, num_bins, num_bins + 1, base=sf))[:-1]
     high = list(sp.logspace(0, num_bins, num_bins + 1, base=sf))[1:]
     #
     # Adding "catch all" bins for anything between 0 - 1 and less than 0
     if self.data_vector[0] < 1.0:
         low.insert(0, 0.0)
         high.insert(0, 1.0)
     if self.data_vector[0] < 0.0:
         low.insert(0, self.data_vector[0])
         high.insert(0, 0.0)
     #
     self.bins = [bin_ for bin_ in zip(low, high)]

    def makeRange(self,start,end,points,log):
        """
        Should be pretty self evident it defines a range:

        - float(start)
        - float(end)
        - int(points)
        - bool(log)
        """
        if log:
            rng = scipy.logspace(scipy.log10(start),scipy.log10(end),points)
        else:
            rng = scipy.linspace(start,end,points)
        return rng

def betweenness_distrib(net, use_weights=True, log=False):
    '''
    Computing the betweenness distribution of a network
    
    Parameters
    ----------
    net : :class:`~nngt.Graph` or subclass
        the network to analyze.
    use_weights : bool, optional (default: True)
        use weighted degrees (do not take the sign into account : all weights
        are positive).
    log : bool, optional (default: False)
        use log-spaced bins.
    
    Returns
    -------
    ncounts : :class:`numpy.array`
        number of nodes in each bin
    nbetw : :class:`numpy.array`
        bins for node betweenness
    ecounts : :class:`numpy.array`
        number of edges in each bin
    ebetw : :class:`numpy.array`
        bins for edge betweenness
    '''
    ia_nbetw, ia_ebetw = net.get_betweenness(use_weights)
    num_nbins, num_ebins = int(len(ia_nbetw) / 50), int(len(ia_ebetw) / 50)
    ra_nbins = sp.linspace(ia_nbetw.min(), ia_nbetw.max(), num_nbins)
    ra_ebins = sp.linspace(ia_ebetw.min(), ia_ebetw.max(), num_ebins)
    if log:
        ra_nbins = sp.logspace(sp.log10(sp.maximum(ia_nbetw.min(),10**-8)),
                               sp.log10(ia_nbetw.max()), num_nbins)
        ra_ebins = sp.logspace(sp.log10(sp.maximum(ia_ebetw.min(),10**-8)),
                               sp.log10(ia_ebetw.max()), num_ebins)
    ncounts,nbetw = sp.histogram(ia_nbetw, ra_nbins)
    ecounts,ebetw = sp.histogram(ia_ebetw, ra_ebins)
    return ncounts, nbetw[:-1], ecounts, ebetw[:-1]