Simulated model output
196105
def output(self):
        """"""
        if self._timeseries is None:
            self.compute()
        output = self._timeseries[:, self.system.output_vars]
        if isinstance(self.system, NetworkModel):
            return self.system._reshape_output(output)
        else:
            return output
32100
def simulate():
    ''''''
    #separate model execution from parameters for safe import from other files
    nest.ResetKernel()

    '''
    Configuration of the simulation kernel by the previously defined time
    resolution used in the simulation. Setting "print_time" to True prints
    the already processed simulation time as well as its percentage of the
    total simulation time.
    '''

    nest.SetKernelStatus({"resolution": dt, "print_time": True,
                          "overwrite_files": True})
    
    print("Building network")
    
    
    '''
    Configuration of the model `iaf_psc_alpha` and `poisson_generator`
    using SetDefaults(). This function expects the model to be the
    inserted as a string and the parameter to be specified in a
    dictionary. All instances of theses models created after this point
    will have the properties specified in the dictionary by default.
    '''
    
    nest.SetDefaults("iaf_psc_alpha", neuron_params)
    nest.SetDefaults("poisson_generator",{"rate": p_rate})
    
    '''
    Creation of the nodes using `Create`. We store the returned handles in
    variables for later reference. Here the excitatory and inhibitory, as
    well as the poisson generator and two spike detectors. The spike
    detectors will later be used to record excitatory and inhibitory
    spikes.
    '''
    
    nodes_ex = nest.Create("iaf_psc_alpha",NE)
    nodes_in = nest.Create("iaf_psc_alpha",NI)
    noise    = nest.Create("poisson_generator")
    espikes  = nest.Create("spike_detector")
    ispikes  = nest.Create("spike_detector")
    
    print("first exc node: {}".format(nodes_ex[0]))
    print("first inh node: {}".format(nodes_in[0]))
    
    '''
    distribute membrane potentials
    '''
    nest.SetStatus(nodes_ex, "V_m",
                   random.rand(len(nodes_ex))*neuron_params["V_th"])
    nest.SetStatus(nodes_in, "V_m",
                   random.rand(len(nodes_in))*neuron_params["V_th"])
    
    '''
    Configuration of the spike detectors recording excitatory and
    inhibitory spikes using `SetStatus`, which expects a list of node
    handles and a list of parameter dictionaries. Setting the variable
    "to_file" to True ensures that the spikes will be recorded in a .gdf
    file starting with the string assigned to label. Setting "withtime"
    and "withgid" to True ensures that each spike is saved to file by
    stating the gid of the spiking neuron and the spike time in one line.
    '''
    
    nest.SetStatus(espikes,[{
                       "label": os.path.join(spike_output_path, label + "-EX"),
                       "withtime": True,
                       "withgid": True,
                       "to_file": True,
                       }])
    
    nest.SetStatus(ispikes,[{
                       "label": os.path.join(spike_output_path, label + "-IN"),
                       "withtime": True,
                       "withgid": True,
                       "to_file": True,}])
    
    print("Connecting devices")
    
    '''
    Definition of a synapse using `CopyModel`, which expects the model
    name of a pre-defined synapse, the name of the customary synapse and
    an optional parameter dictionary. The parameters defined in the
    dictionary will be the default parameter for the customary
    synapse. Here we define one synapse for the excitatory and one for the
    inhibitory connections giving the previously defined weights and equal
    delays.
    '''
    
    nest.CopyModel("static_synapse","excitatory",{"weight":J_ex, "delay":delay})
    nest.CopyModel("static_synapse","inhibitory",{"weight":J_in, "delay":delay})
    
    '''
    Connecting the previously defined poisson generator to the excitatory
    and inhibitory neurons using the excitatory synapse. Since the poisson
    generator is connected to all neurons in the population the default
    rule ('all_to_all') of Connect() is used. The synaptic properties are
    inserted via syn_spec which expects a dictionary when defining
    multiple variables or a string when simply using a pre-defined
    synapse.
    '''
    
    if Poisson:
      nest.Connect(noise,nodes_ex, 'all_to_all', "excitatory")
      nest.Connect(noise,nodes_in,'all_to_all', "excitatory")
    
    '''
    Connecting the first N_neurons nodes of the excitatory and inhibitory
    population to the associated spike detectors using excitatory
    synapses. Here the same shortcut for the specification of the synapse
    as defined above is used.
    '''
    
    nest.Connect(nodes_ex,espikes, 'all_to_all', "excitatory")
    nest.Connect(nodes_in,ispikes, 'all_to_all', "excitatory")
    
    print("Connecting network")
    
    print("Excitatory connections")
    
    '''
    Connecting the excitatory population to all neurons using the
    pre-defined excitatory synapse. Beforehand, the connection parameter
    are defined in a dictionary. Here we use the connection rule
    'fixed_indegree', which requires the definition of the indegree. Since
    the synapse specification is reduced to assigning the pre-defined
    excitatory synapse it suffices to insert a string.
    '''
    
    conn_params_ex = {'rule': 'fixed_indegree', 'indegree': CE}
    nest.Connect(nodes_ex, nodes_ex+nodes_in, conn_params_ex, "excitatory")
    
    print("Inhibitory connections")
    
    '''
    Connecting the inhibitory population to all neurons using the
    pre-defined inhibitory synapse. The connection parameter as well as
    the synapse paramtere are defined analogously to the connection from
    the excitatory population defined above.
    '''
    
    conn_params_in = {'rule': 'fixed_indegree', 'indegree': CI}
    nest.Connect(nodes_in, nodes_ex+nodes_in, conn_params_in, "inhibitory")
    
    
    '''
    Storage of the time point after the buildup of the network in a
    variable.
    '''
    
    endbuild=time.time()
    
    '''
    Simulation of the network.
    '''
    
    print("Simulating")
    
    nest.Simulate(simtime)
    
    '''
    Storage of the time point after the simulation of the network in a
    variable.
    '''
    
    endsimulate= time.time()
    
    '''
    Reading out the total number of spikes received from the spike
    detector connected to the excitatory population and the inhibitory
    population.
    '''
    
    events_ex = nest.GetStatus(espikes,"n_events")[0]
    events_in = nest.GetStatus(ispikes,"n_events")[0]
    
    '''
    Calculation of the average firing rate of the excitatory and the
    inhibitory neurons by dividing the total number of recorded spikes by
    the number of neurons recorded from and the simulation time. The
    multiplication by 1000.0 converts the unit 1/ms to 1/s=Hz.
    '''
    
    rate_ex   = events_ex/simtime*1000.0/N_neurons
    rate_in   = events_in/simtime*1000.0/N_neurons
    
    '''
    Reading out the number of connections established using the excitatory
    and inhibitory synapse model. The numbers are summed up resulting in
    the total number of synapses.
    '''
    
    num_synapses = nest.GetDefaults("excitatory")["num_connections"]+\
    nest.GetDefaults("inhibitory")["num_connections"]
    
    '''
    Establishing the time it took to build and simulate the network by
    taking the difference of the pre-defined time variables.
    '''
    
    build_time = endbuild-startbuild
    sim_time   = endsimulate-endbuild
    
    '''
    Printing the network properties, firing rates and building times.
    '''
    
    print("Brunel network simulation (Python)")
    print("Number of neurons : {0}".format(N_neurons))
    print("Number of synapses: {0}".format(num_synapses))
    print("       Exitatory  : {0}".format(int(CE * N_neurons) + N_neurons))
    print("       Inhibitory : {0}".format(int(CI * N_neurons)))
    print("Excitatory rate   : %.2f Hz" % rate_ex)
    print("Inhibitory rate   : %.2f Hz" % rate_in)
    print("Building time     : %.2f s" % build_time)
    print("Simulation time   : %.2f s" % sim_time)
    
    '''
    Plot a raster of the excitatory neurons and a histogram.
    '''
    
    if False:
        nest.raster_plot.from_device(espikes, hist=True)
        nest.raster_plot.from_device(ispikes, hist=True)
        nest.raster_plot.show()
285064
def computeOutput(self):

    """"""

    # todo: This operation can be sped up by:
    #  1.)  Pre-allocating space for the currentOutput
    #  2.)  Making predictedState and activeState of type 'float32' up front
    #  3.)  Using logical_or(self.predictedState['t'], self.activeState['t'],
    #          self.currentOutput)
    #
    self.currentOutput = \
        logical_or(self.predictedState['t'], self.activeState['t'])

    return self.currentOutput.reshape(-1).astype('float32')
152176
def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False):
    """"""
    Q_signal = 4
    import GPy
    import numpy as np
    np.random.seed(3000)

    k = GPy.kern.Matern32(Q_signal, 1., lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1)
    for i in range(Q_signal):
        k += GPy.kern.PeriodicExponential(1, variance=1., active_dims=[i], period=3., lower=-2, upper=6)
    t = np.c_[[np.linspace(-1, 5, N) for _ in range(Q_signal)]].T
    K = k.K(t)
    s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:, :, None]

    Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS)

    slist = [sS, s1, s2, s3]
    slist_names = ["sS", "s1", "s2", "s3"]
    Ylist = [Y1, Y2, Y3]

    if plot_sim:
        from matplotlib import pyplot as plt
        import matplotlib.cm as cm
        import itertools
        fig = plt.figure("MRD Simulation Data", figsize=(8, 6))
        fig.clf()
        ax = fig.add_subplot(2, 1, 1)
        labls = slist_names
        for S, lab in zip(slist, labls):
            ax.plot(S, label=lab)
        ax.legend()
        for i, Y in enumerate(Ylist):
            ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
            ax.imshow(Y, aspect='auto', cmap=cm.gray)  # @UndefinedVariable
            ax.set_title("Y{}".format(i + 1))
        plt.draw()
        plt.tight_layout()

    return slist, [S1, S2, S3], Ylist
1199452
def generate(self):
        '''
        
        '''
        observed_arr = None
        for result_tuple in self.__feature_generator.generate():
            observed_arr = result_tuple[0]
            break

        if self.noise_sampler is not None:
            self.noise_sampler.output_shape = observed_arr.shape
            observed_arr += self.noise_sampler.generate()

        observed_arr = observed_arr.astype(float)
        if self.__norm_mode == "z_score":
            if observed_arr.std() != 0:
                observed_arr = (observed_arr - observed_arr.mean()) / observed_arr.std()
        elif self.__norm_mode == "min_max":
            if (observed_arr.max() - observed_arr.min()) != 0:
                observed_arr = (observed_arr - observed_arr.min()) / (observed_arr.max() - observed_arr.min())
        elif self.__norm_mode == "tanh":
            observed_arr = np.tanh(observed_arr)
        
        return observed_arr
558780
def _output(cls, fluents: Sequence[FluentPair]) -> Sequence[tf.Tensor]:
        ''''''
        output = []
        for _, fluent in fluents:
            tensor = fluent.tensor
            if tensor.dtype != tf.float32:
                tensor = tf.cast(tensor, tf.float32)
            output.append(tensor)
        return tuple(output)
556507
def mu(self):
        """"""
        mu = self._models[0].mu
        assert all([mu == model.mu for model in self._models])
        return mu
856904
def estimate(self, data, full_output=False, **kwargs):
        

        # Number of model comparisons can be specified in the configuration.
        num_model_comparisons = self._configuration.get("estimate", {}).get(
            "num_model_comparisons", self.grid_points.size)
        # If it's a fraction, we need to convert that to an integer.
        if 1 > num_model_comparisons > 0:
            num_model_comparisons *= self.grid_points.size

        # If the num_model_comparison is provided as a keyword argument, use it.
        num_model_comparisons = kwargs.pop("num_model_comparisons",
            int(num_model_comparisons))

        logger.debug("Number of model comparisons to make for initial estimate:"
            " {0}".format(num_model_comparisons))
        
        # Match the data to the model channels.
        matched_channels, missing_channels, ignore_parameters \
            = self._match_channels_to_data(data)

        logger.debug("Matched channels: {0}, missing channels: {1}, ignore "
            "parameters: {2}".format(matched_channels, missing_channels,
                ignore_parameters))

        # Load the intensities
        t = time()
        s = self.grid_points.size/num_model_comparisons # step size
        grid_points = self.grid_points[::s]
        intensities = np.memmap(
            self._configuration["model_grid"]["intensities"], dtype="float32",
            mode="r", shape=(self.grid_points.size, self.wavelengths.size))[::s]
        logger.debug("Took {:.0f} seconds to load and slice intensities".format(
            time() - t))
        # Which matched, data channel has the highest S/N?
        # (This channel will be used to estimate astrophysical parameters)
        data, pixels_affected = self._apply_data_mask(data)
        median_snr = dict(zip(matched_channels,
            [np.nanmedian(spec.flux/(spec.variance**0.5)) for spec in data]))
        median_snr.pop(None, None) # Remove unmatched data spectra

        ccf_channel = self._configuration.get("settings", {}).get("ccf_channel",
            max(median_snr, key=median_snr.get))
        if ccf_channel not in matched_channels:
            logger.warn("Ignoring CCF channel {0} because it was not a matched"
                " channel".format(ccf_channel))
            ccf_channel = max(median_snr, key=median_snr.get)

        logger.debug("Channel with peak SNR is {0}".format(ccf_channel))

        # Are there *any* continuum parameters in any matched channel?
        any_continuum_parameters = any(map(lambda s: s.startswith("continuum_"),
            set(self.parameters).difference(ignore_parameters)))

        # [TODO]: CCF MASK
        # [TODO]: Don't require CCF if we have only continuum parameters.

        z_limits = self._configuration["settings"].get("ccf_z_limits", None)

        theta = {} # Dictionary for the estimated model parameters.
        best_grid_index = None
        c = speed_of_light.to("km/s").value
        for matched_channel, spectrum in zip(matched_channels, data):
            if matched_channel is None: continue

            # Do we need todo cross-correlation for this channel?
            # We do if there are redshift parameters for this channel,
            # or if there is a global redshift or global continuum parameters
            #   and this channel is the highest S/N.
            if "z_{}".format(matched_channel) in self.parameters \
            or ((any_continuum_parameters or "z" in self.parameters) \
            and matched_channel == ccf_channel):

                # Get the continuum degree for this channel.
                continuum_degree = self._configuration["model"].get("continuum",
                    { matched_channel: -1 })[matched_channel]

                logger.debug("Perfoming CCF on {0} channel with a continuum "
                    "degree of {1}".format(matched_channel, continuum_degree))

                # Get model wavelength indices that match the data.
                # get the points that are in the mask, and within the spectrum
                # limits

                # TODO: Make this CCF not model mask.
                idx = np.where(self._model_mask() \
                    * (self.wavelengths >= spectrum.disp[0]) \
                    * (spectrum.disp[-1] >= self.wavelengths))[0]

                v, v_err, R = spectrum.cross_correlate(
                    (self.wavelengths[idx], intensities[:, idx]),
                    #(self.wavelengths, intensities),
                    continuum_degree=continuum_degree, z_limits=z_limits)

                # Identify the best point by the CCF peak.
                best = np.nanargmax(R)

                # Now, why did we do CCF in this channel? Which model parameters
                # should be updated?
                if "z_{}".format(matched_channel) in self.parameters:
                    theta["z_{}".format(matched_channel)] = v[best] / c
                elif "z" in self.parameters: 
                    # If there is a global redshift, update it.
                    theta["z"] = v[best] / c

                    # Continuum parameters will be updated later, so that each
                    # channel is checked to see if it has the highest S/N,
                    # otherwise we might be trying to calculate continuum
                    # parameters when we haven't done CCF on the highest S/N
                    # spectra yet.

                if matched_channel == ccf_channel:
                    # Update astrophysical parameters.
                    theta.update(dict(zip(grid_points.dtype.names,
                        grid_points[best])))
                    best_grid_index = best

        # If there are continuum parameters, calculate them from the best point.
        if any_continuum_parameters:
            for matched_channel, spectrum in zip(matched_channels, data):
                if matched_channel is None: continue

                # The template spectra at the best point needs to be
                # redshifted to the data, and then continuum coefficients
                # calculated from that.

                # Get the continuum degree for this channel.
                continuum_degree = self._configuration["model"].get("continuum",
                    { matched_channel: -1 })[matched_channel]

            
                # Get model wavelength indices that match the data.
                idx = np.clip(self.wavelengths.searchsorted(
                    [spectrum.disp[0], spectrum.disp[-1]]) + [0, 1],
                    0, self.wavelengths.size)

                # Redshift and bin the spectrum.
                z = theta.get("z_{}".format(matched_channel), theta.get("z", 0))

                best_intensities \
                    = np.copy(intensities[best_grid_index, idx[0]:idx[1]]).flatten()

                # Apply model mask.
                model_mask = self._model_mask(self.wavelengths[idx[0]:idx[1]])
                best_intensities[~model_mask] = np.nan

                best_intensities = best_intensities * specutils.sample.resample(
                    self.wavelengths[idx[0]:idx[1]] * (1 + z), spectrum.disp)
                
                # Calculate the continuum coefficients for this channel.
                continuum = spectrum.flux/best_intensities
                finite = np.isfinite(continuum)

                try:
                    coefficients = np.polyfit(
                        spectrum.disp[finite], continuum[finite], continuum_degree,
                        )#w=spectrum.ivariance[finite])
                except np.linalg.linalg.LinAlgError:
                    logger.exception("Exception in initial polynomial fit")
                    coefficients = np.polyfit(spectrum.disp[finite], continuum[finite],
                        continuum_degree)

                # They go into theta backwards. such that coefficients[-1] is
                # continuum_{name}_0
                theta.update(dict(zip(
                    ["continuum_{0}_{1}".format(matched_channel, i) \
                        for i in range(continuum_degree + 1)],
                    coefficients[::-1]
                )))

        # Remaining parameters could be: resolving power, outlier pixels,
        # underestimated variance.
        remaining_parameters = set(self.parameters)\
            .difference(ignore_parameters)\
            .difference(theta)

        if remaining_parameters:
            logger.debug("Remaining parameters to estimate: {0}. For these we "
                "will just assume reasonable initial values.".format(
                remaining_parameters))

            for parameter in remaining_parameters:
                if parameter == "resolution" \
                or parameter.startswith("resolution_"):

                    if parameter.startswith("resolution_"):
                        spectra = [data[matched_channels.index(
                            parameter.split("_")[1])]]
                    else:
                        spectra = [s for s in data if s is not None]

                    R = [s.disp.mean()/np.diff(s.disp).mean() for s in spectra]

                    # Assume oversampling rate of ~5.
                    theta.update({ parameter: np.median(R)/5.})

                elif parameter == "ln_f" or parameter.startswith("ln_f_"):
                    theta.update({ parameter: 0.5 }) # Not overestimated.

                elif parameter in ("Po", "Vo"):
                    theta.update({
                        "Po": 0.01, # 1% outlier pixels.
                        "Vo": np.mean([np.nanmedian(s.variance) for s in data]),
                    })

        logger.info("Initial estimate: {}".format(theta))
        # Having full_output = True means return the best spectra estimate.
        if full_output:

            # Create model fluxes and calculate some metric.
            __intensities = np.copy(intensities[best_grid_index])

            # Apply model masks.
            __intensities[~self._model_mask()] = np.nan

            chi_sq, dof, model_fluxes = self._chi_sq(theta, data,
                __intensities=__intensities, __no_precomputed_binning=True)
            del intensities

            return (theta, chi_sq, dof, model_fluxes)

        # Delete the reference to intensities
        del intensities
        return theta
496616
def main():
    
    input_tif = "../tests/data/Jamaica_dem.tif"
    output_tif = "../tests/data/tmp_results/log_dem.tif"
    rst = RasterUtilClass.read_raster(input_tif)
    # raster data (with noDataValue as numpy.nan) as numpy array
    rst_valid = rst.validValues
    output_data = np.log(rst_valid)
    # write output raster
    RasterUtilClass.write_gtiff_file(output_tif, rst.nRows, rst.nCols, output_data, rst.geotrans,
                                     rst.srs, rst.noDataValue, rst.dataType)
298380
def estimate(args, **kwargs):
    

    expected_output_files = kwargs.pop("expected_output_files", None)
    if not expected_output_files:
        expected_output_files = ["estimate.pkl"]
        if args.plotting:
            expected_output_files.extend(
                ["projection-estimate.{}".format(args.plot_format)])
        
    model, data, metadata = _pre_solving(args, expected_output_files)

    try:
        theta, chisq, dof, model_fluxes = model.estimate(data, full_output=True,
            debug=args.debug)

    except:
        logger.exception("Failed to estimate model parameters")
        raise

    logger.info("Estimated model parameters are:")
    _announce_theta(theta)
    logger.info("With a chi-sq value of {0:.1f} (reduced {1:.1f}; DOF {2:.1f})"\
        .format(chisq, chisq/dof, dof))

    metadata["estimated"] = {
        "theta": theta,
        "chi_sq": chisq,
        "dof": dof,
        "r_chi_sq": chisq/dof
    }

    if args.plotting:
        fig = sick.plot.spectrum(data, model_flux=model_fluxes)
        filename = _prefix(args, "projection-estimate.{}".format(
            args.plot_format))
        fig.savefig(filename)
        logger.info("Created figure {}".format(filename))

    if kwargs.pop("__return_result", False):
        return (model, data, metadata, theta)

    # Write the result to file.
    _write_output(_prefix(args, "estimate.pkl"), metadata)
    return None
1209718
private double simulateParameters(double[] curparams) {
    // Initialize alpha, beta
    for(int i = 0; i < numfit; i++) {
      Arrays.fill(alpha[i], 0.);
    }
    Arrays.fill(beta, 0.);

    double newchisq = 0.0;

    // Simulation loop over all data
    for(int di = 0; di < x.length; di++) {
      FittingFunctionResult res = func.eval(x[di], curparams);
      // compute inverse squared standard deviation of the point (confidence?)
      double sigma2inv = 1.0 / (s[di] * s[di]);
      double deltay = y[di] - res.y;
      // i2 and j2 are the indices that only count the params with dofit true!
      for(int i = 0, i2 = 0; i < numfit; i++) {
        if(dofit[i]) {
          double wt = res.gradients[i] * sigma2inv;
          // fill only half of the matrix, use symmetry below to complete the
          // remainder.
          for(int j = 0, j2 = 0; j <= i; j++) {
            if(dofit[j]) {
              alpha[i2][j2++] += wt * res.gradients[j];
            }
          }
          beta[i2++] += deltay * wt;
        }
      }
      newchisq += deltay * deltay * sigma2inv;
    }
    // fill symmetric side of matrix
    for(int i = 1; i < numfit; i++) {
      for(int j = i + 1; j < numfit; j++) {
        alpha[i][j] = alpha[j][i];
      }
    }

    return newchisq;
  }
128327
def _output_cdf_plot(core_result, spid, models, options, fit_results):
    """"""

    # CDF
    x = core_result['y'].values
    df = emp.empirical_cdf(x)
    df.columns = ['x', 'empirical']

    def calc_func(model, df, shapes):
        return eval("mod.%s.cdf(df['x'], *shapes)" % model)

    plot_exec_str = "ax.step(df['x'], emp, color='k', lw=3);ax.set_ylim(top=1)"

    _save_table_and_plot(spid, models, options, fit_results, 'data_pred_cdf',
                         df, calc_func, plot_exec_str)
844664
def _xfs_estimate_output(out):
    '''
    
    '''
    spc = re.compile(r"\s+")
    data = {}
    for line in [l for l in out.split("\n") if l.strip()][1:]:
        directory, bsize, blocks, megabytes, logsize = spc.sub(" ", line).split(" ")
        data[directory] = {
            'block _size': bsize,
            'blocks': blocks,
            'megabytes': megabytes,
            'logsize': logsize,
        }

    return data
590900
def writeOutput(self, filename, samples, srcFs, targetFs):
        
        import librosa
        inputLength = samples.shape[-1]
        if srcFs != targetFs:
            if inputLength < targetFs:
                logger.info("Input signal is too short (" + str(inputLength) +
                            " samples) for resampling to " + str(targetFs) + "Hz")
                outputSamples = samples
                targetFs = srcFs
            else:
                logger.info("Resampling " + str(inputLength) + " samples from " + str(srcFs) + "Hz to " +
                            str(targetFs) + "Hz")
                outputSamples = librosa.resample(samples, srcFs, targetFs, res_type='kaiser_fast')
        else:
            outputSamples = samples
        logger.info("Writing output to " + filename)
        maxv = np.iinfo(np.int32).max
        librosa.output.write_wav(filename, (outputSamples * maxv).astype(np.int32), targetFs)
        logger.info("Output written to " + filename)
958560
def simulate(self, partition, outdir, jobhandler=default_jobhandler, batchsize=1, **kwargs):
        

        results = self.get_partition_results(partition)
        DEFAULT_DNA_MODEL = 'GTR'
        DEFAULT_PROTEIN_MODEL = 'LG08'

        # Collect argument list
        args = [None] * len(self.collection)

        for result in results:
            if len(result['partitions']) > 1:
                places = dict((j,i) for (i,j) in enumerate(rec.name for rec in self.collection.records))
                for partition in result['partitions'].values():
                    place = places[partition['name']]
                    model = partition.get('model')
                    freqs = partition.get('frequencies')
                    rates = partition.get('rates')
                    alpha = partition.get('alpha')
                    tree  = str(result['ml_tree'])
                    if model is None:
                        model = DEFAULT_DNA_MODEL if self.collection[place].is_dna() else DEFAULT_PROTEIN_MODEL
                    if freqs is not None:
                        freqs = smooth_freqs(freqs)
                    args[place] = (len(self.collection[place]),
                                   model_translate(model),
                                   freqs,
                                   alpha,
                                   tree,
                                   rates)
            else:
                model = result['partitions']['0'].get('model')
                freqs = result['partitions']['0'].get('frequencies')
                rates = result['partitions']['0'].get('rates')
                alpha = result['partitions']['0'].get('alpha')
                tree  = str(result['ml_tree'])
                if freqs is not None:
                    freqs = smooth_freqs(freqs)
                use_default_model = (model is None)
                for i in range(len(self.collection)):
                    if use_default_model:
                        model = DEFAULT_DNA_MODEL if self.collection[i].is_dna() else DEFAULT_PROTEIN_MODEL
                    args[i] = (len(self.collection[i]),
                               model_translate(model),
                               freqs,
                               alpha,
                               tree,
                               rates)

        # Distribute work
        msg = 'Simulating'
        map_result = jobhandler(tasks.simulate_task, args, msg, batchsize)

        # Process results
        for i, result in enumerate(map_result):
            orig = self.collection[i]
            simseqs = gapmask(result.items(), orig.get_sequences())
            al = Alignment(simseqs, alphabet=('protein' if orig.is_protein() else 'dna'))
            outfile = os.path.join(outdir, orig.name + '.phy')
            al.write_alignment(outfile, 'phylip', True)
735952
def draw(self):
        '''
        
        '''
        observed_arr = self.extract_conditions()
        conv_arr = self.inference(observed_arr)

        if self.__conditon_noise_sampler is not None:
            self.__conditon_noise_sampler.output_shape = conv_arr.shape
            noise_arr = self.__conditon_noise_sampler.generate()
            conv_arr += noise_arr

        deconv_arr = self.__deconvolution_model.inference(conv_arr)
        return np.concatenate((deconv_arr, observed_arr), axis=1)
730870
def _prepare_output_multi(self, model):
        """"""
        model_name = model.__name__
        current_path = os.path.join(self._output_path, '{model}.{extension}'.format(
            model=model_name,
            extension=self.EXTENSION,
        ))
        self._outfile = codecs.open(current_path, 'w', encoding='utf-8')
        print('Dumping {model} to {file}'.format(model=model_name, file=current_path))